Differential D58 Diff 240 src/EvtGenExternal/EvtPhotosEngine.cpp

Changeset 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