Differential D58 Diff 252 src/EvtGenExternal/EvtTauolaEngine.cpp

Changeset View

Standalone View

src/EvtGenExternal/EvtTauolaEngine.cpp

Show All 37 Lines
#include <sstream>		#include <sstream>
#include <string>		#include <string>

using std::endl;		using std::endl;

EvtTauolaEngine::EvtTauolaEngine( bool useEvtGenRandom )		EvtTauolaEngine::EvtTauolaEngine( bool useEvtGenRandom )
{		{
// PDG standard code integer ID for tau particle		// PDG standard code integer ID for tau particle
_tauPDG = 15;		m_tauPDG = 15;
// Number of possible decay modes in Tauola		// Number of possible decay modes in Tauola
_nTauolaModes = 22;		m_nTauolaModes = 22;

EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Setting up TAUOLA." << endl;		EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Setting up TAUOLA." << endl;

// These three lines are not really necessary since they are the default.		// These three lines are not really necessary since they are the default.
// But they are here so that we know what the initial conditions are.		// But they are here so that we know what the initial conditions are.
Tauolapp::Tauola::setDecayingParticle( _tauPDG ); // tau PDG code		Tauolapp::Tauola::setDecayingParticle( m_tauPDG ); // tau PDG code
Tauolapp::Tauola::setSameParticleDecayMode(		Tauolapp::Tauola::setSameParticleDecayMode(
Tauolapp::Tauola::All ); // all modes allowed		Tauolapp::Tauola::All ); // all modes allowed
Tauolapp::Tauola::setOppositeParticleDecayMode(		Tauolapp::Tauola::setOppositeParticleDecayMode(
Tauolapp::Tauola::All ); // all modes allowed		Tauolapp::Tauola::All ); // all modes allowed

// Limit the number of warnings printed out. Can't choose zero here, unfortunately		// Limit the number of warnings printed out. Can't choose zero here, unfortunately
Tauolapp::Log::SetWarningLimit( 1 );		Tauolapp::Log::SetWarningLimit( 1 );

// Initial the Tauola external generator		// Initial the Tauola external generator
if ( useEvtGenRandom == true ) {		if ( useEvtGenRandom == true ) {
EvtGenReport( EVTGEN_INFO, "EvtGen" )		EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "Using EvtGen random number engine also for Tauola++" << endl;		<< "Using EvtGen random number engine also for Tauola++" << endl;

Tauolapp::Tauola::setRandomGenerator( EvtRandom::Flat );		Tauolapp::Tauola::setRandomGenerator( EvtRandom::Flat );
}		}

// Use the BaBar-tuned chiral current calculations by default. Can be changed using the		// Use the BaBar-tuned chiral current calculations by default. Can be changed using the
// TauolaCurrentOption keyword in decay files		// TauolaCurrentOption keyword in decay files
Tauolapp::Tauola::setNewCurrents( 1 );		Tauolapp::Tauola::setNewCurrents( 1 );

Tauolapp::Tauola::initialize();		Tauolapp::Tauola::initialize();

// Initialise various default parameters		// Initialise various default parameters
// Neutral and charged spin propagator choices		// Neutral and charged spin propagator choices
_neutPropType = 0;		m_neutPropType = 0;
_posPropType = 0;		m_posPropType = 0;
_negPropType = 0;		m_negPropType = 0;

// Set-up possible decay modes _after_ we have read the (user) decay file		// Set-up possible decay modes _after_ we have read the (user) decay file
_initialised = false;		m_initialised = false;
}		}

void EvtTauolaEngine::initialise()		void EvtTauolaEngine::initialise()
{		{
// Set up all possible tau decay modes.		// Set up all possible tau decay modes.
// This should be done just before the first doDecay() call,		// This should be done just before the first doDecay() call,
// since we want to make sure that any decay.dec files are processed		// since we want to make sure that any decay.dec files are processed
// first to get lists of particle modes and their alias definitions		// first to get lists of particle modes and their alias definitions
// (for creating EvtParticles with the right history information).		// (for creating EvtParticles with the right history information).

if ( _initialised == false ) {		if ( m_initialised == false ) {
this->setUpPossibleTauModes();		this->setUpPossibleTauModes();
this->setOtherParameters();		this->setOtherParameters();

_initialised = true;		m_initialised = true;
}		}
}		}

void EvtTauolaEngine::setUpPossibleTauModes()		void EvtTauolaEngine::setUpPossibleTauModes()
{		{
// Get the decay table list defined by the decay.dec files.		// Get the decay table list defined by the decay.dec files.
// Only look for the first tau particle decay mode definitions with the Tauola name,		// Only look for the first tau particle decay mode definitions with the Tauola name,
// since that generator only allows the same BFs for both tau+ and tau- decays.		// since that generator only allows the same BFs for both tau+ and tau- decays.
// We can not choose a specific tau decay event-by-event, since this is		// We can not choose a specific tau decay event-by-event, since this is
// only possible before we call Tauola::initialize().		// only possible before we call Tauola::initialize().
// Otherwise, we could have selected a random mode ourselves for tau- and tau+		// Otherwise, we could have selected a random mode ourselves for tau- and tau+
// separately (via selecting a random number and comparing it to be less than		// separately (via selecting a random number and comparing it to be less than
// the cumulative BF) for each event.		// the cumulative BF) for each event.

int nPDL = EvtPDL::entries();		int nPDL = EvtPDL::entries();
int iPDL( 0 );		int iPDL( 0 );

bool gotAnyTauolaModes( false );		bool gotAnyTauolaModes( false );

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

if ( abs( PDGId ) == _tauPDG && gotAnyTauolaModes == false ) {		if ( abs( PDGId ) == m_tauPDG && gotAnyTauolaModes == false ) {
int aliasInt = particleId.getAlias();		int aliasInt = particleId.getAlias();

// Get the list of decay modes for this tau particle (alias)		// Get the list of decay modes for this tau particle (alias)
int nModes = EvtDecayTable::getInstance()->getNModes( aliasInt );		int nModes = EvtDecayTable::getInstance()->getNModes( aliasInt );
int iMode( 0 ), iTauMode( 0 );		int iMode( 0 ), iTauMode( 0 );

// Vector to store tau mode branching fractions.		// Vector to store tau mode branching fractions.
// The size of this vector equals the total number of possible		// The size of this vector equals the total number of possible
// Tauola decay modes. Initialise all BFs to zero.		// Tauola decay modes. Initialise all BFs to zero.
std::vector<double> tauolaModeBFs( _nTauolaModes );		std::vector<double> tauolaModeBFs( m_nTauolaModes );

for ( iTauMode = 0; iTauMode < _nTauolaModes; iTauMode++ ) {		for ( iTauMode = 0; iTauMode < m_nTauolaModes; iTauMode++ ) {
tauolaModeBFs[iTauMode] = 0.0;		tauolaModeBFs[iTauMode] = 0.0;
}		}

double totalTauModeBF( 0.0 );		double totalTauModeBF( 0.0 );

int nNonTauolaModes( 0 );		int nNonTauolaModes( 0 );

// Loop through each decay mode		// Loop through each decay mode
for ( iMode = 0; iMode < nModes; iMode++ ) {		for ( iMode = 0; iMode < nModes; iMode++ ) {
EvtDecayBase* decayModel =		EvtDecayBase* decayModel =
EvtDecayTable::getInstance()->findDecayModel( aliasInt,		EvtDecayTable::getInstance()->findDecayModel( aliasInt,
iMode );		iMode );
if ( decayModel ) {		if ( decayModel ) {
// Check that the decay model name matches TAUOLA		// Check that the decay model name matches TAUOLA
std::string modelName = decayModel->getName();		std::string modelName = decayModel->getName();
if ( modelName == "TAUOLA" ) {		if ( modelName == "TAUOLA" ) {
if ( gotAnyTauolaModes == false ) {		if ( gotAnyTauolaModes == false ) {
gotAnyTauolaModes = true;		gotAnyTauolaModes = true;
}		}

// Extract the decay mode integer type and branching fraction		// Extract the decay mode integer type and branching fraction
double BF = decayModel->getBranchingFraction();		double BF = decayModel->getBranchingFraction();
int modeArrayInt = this->getModeInt( decayModel ) - 1;		int modeArrayInt = this->getModeInt( decayModel ) - 1;

if ( modeArrayInt >= 0 && modeArrayInt < _nTauolaModes ) {		if ( modeArrayInt >= 0 && modeArrayInt < m_nTauolaModes ) {
tauolaModeBFs[modeArrayInt] = BF;		tauolaModeBFs[modeArrayInt] = BF;
totalTauModeBF += BF;		totalTauModeBF += BF;
}		}

} else {		} else {
nNonTauolaModes++;		nNonTauolaModes++;
}		}

Show All 10 Lines	for ( iPDL = 0; iPDL < nPDL; iPDL++ ) {

// Normalise all (non-zero) tau mode BFs to sum up to 1.0, and		// Normalise all (non-zero) tau mode BFs to sum up to 1.0, and
// let Tauola know about these normalised branching fractions		// let Tauola know about these normalised branching fractions
if ( totalTauModeBF > 0.0 ) {		if ( totalTauModeBF > 0.0 ) {
EvtGenReport( EVTGEN_INFO, "EvtGen" )		EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "Setting TAUOLA BF modes using the definitions for the particle "		<< "Setting TAUOLA BF modes using the definitions for the particle "
<< EvtPDL::name( particleId ) << endl;		<< EvtPDL::name( particleId ) << endl;

for ( iTauMode = 0; iTauMode < _nTauolaModes; iTauMode++ ) {		for ( iTauMode = 0; iTauMode < m_nTauolaModes; iTauMode++ ) {
tauolaModeBFs[iTauMode] /= totalTauModeBF;		tauolaModeBFs[iTauMode] /= totalTauModeBF;
double modeBF = tauolaModeBFs[iTauMode];		double modeBF = tauolaModeBFs[iTauMode];
EvtGenReport( EVTGEN_INFO, "EvtGen" )		EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "Setting TAUOLA BF for mode " << iTauMode + 1		<< "Setting TAUOLA BF for mode " << iTauMode + 1
<< " = " << modeBF << endl;		<< " = " << modeBF << endl;
Tauolapp::Tauola::setTauBr( iTauMode + 1, modeBF );		Tauolapp::Tauola::setTauBr( iTauMode + 1, modeBF );
}		}

Show All 27 Lines	void EvtTauolaEngine::setOtherParameters()
// Set other Tauola parameters using the "Defined" keyword in the decay file. If any of		// Set other Tauola parameters using the "Defined" keyword in the decay file. If any of
// these are not found in the decay file, then default values are assumed/kept		// these are not found in the decay file, then default values are assumed/kept

// 1) TauolaNeutralProp: Specify the neutral propagator type used for spin matrix calculations		// 1) TauolaNeutralProp: Specify the neutral propagator type used for spin matrix calculations
// "Z" (default), "Gamma", "Higgs" (H0), "PseudoHiggs" (A0), "MixedHiggs" (A0/H0)		// "Z" (default), "Gamma", "Higgs" (H0), "PseudoHiggs" (A0), "MixedHiggs" (A0/H0)
int iErr( 0 );		int iErr( 0 );
std::string neutPropName = EvtSymTable::get( "TauolaNeutralProp", iErr );		std::string neutPropName = EvtSymTable::get( "TauolaNeutralProp", iErr );
if ( neutPropName == "Z0" \|\| neutPropName == "Z" ) {		if ( neutPropName == "Z0" \|\| neutPropName == "Z" ) {
_neutPropType = Tauolapp::TauolaParticle::Z0;		m_neutPropType = Tauolapp::TauolaParticle::Z0;
} else if ( neutPropName == "Gamma" ) {		} else if ( neutPropName == "Gamma" ) {
_neutPropType = Tauolapp::TauolaParticle::GAMMA;		m_neutPropType = Tauolapp::TauolaParticle::GAMMA;
} else if ( neutPropName == "Higgs" ) {		} else if ( neutPropName == "Higgs" ) {
_neutPropType = Tauolapp::TauolaParticle::HIGGS;		m_neutPropType = Tauolapp::TauolaParticle::HIGGS;
} else if ( neutPropName == "PseudoHiggs" ) {		} else if ( neutPropName == "PseudoHiggs" ) {
_neutPropType = Tauolapp::TauolaParticle::HIGGS_A;		m_neutPropType = Tauolapp::TauolaParticle::HIGGS_A;
} else if ( neutPropName == "MixedHiggs" ) {		} else if ( neutPropName == "MixedHiggs" ) {
_neutPropType = Tauolapp::Tauola::getHiggsScalarPseudoscalarPDG();		m_neutPropType = Tauolapp::Tauola::getHiggsScalarPseudoscalarPDG();
}		}

if ( _neutPropType != 0 ) {		if ( m_neutPropType != 0 ) {
EvtGenReport( EVTGEN_INFO, "EvtGen" )		EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "TAUOLA neutral spin propagator PDG id set to " << _neutPropType		<< "TAUOLA neutral spin propagator PDG id set to " << m_neutPropType
<< endl;		<< endl;
}		}

// 2) TauolaChargedProp: Specify the charged propagator type used for spin matrix calculations		// 2) TauolaChargedProp: Specify the charged propagator type used for spin matrix calculations
// "W" (default), "Higgs" (H+/H-)		// "W" (default), "Higgs" (H+/H-)
std::string chargedPropName = EvtSymTable::get( "TauolaChargedProp", iErr );		std::string chargedPropName = EvtSymTable::get( "TauolaChargedProp", iErr );
if ( chargedPropName == "W" ) {		if ( chargedPropName == "W" ) {
_negPropType = Tauolapp::TauolaParticle::W_MINUS;		m_negPropType = Tauolapp::TauolaParticle::W_MINUS;
_posPropType = Tauolapp::TauolaParticle::W_PLUS;		m_posPropType = Tauolapp::TauolaParticle::W_PLUS;
} else if ( chargedPropName == "Higgs" ) {		} else if ( chargedPropName == "Higgs" ) {
_negPropType = Tauolapp::TauolaParticle::HIGGS_MINUS;		m_negPropType = Tauolapp::TauolaParticle::HIGGS_MINUS;
_posPropType = Tauolapp::TauolaParticle::HIGGS_PLUS;		m_posPropType = Tauolapp::TauolaParticle::HIGGS_PLUS;
}		}

if ( _negPropType != 0 ) {		if ( m_negPropType != 0 ) {
EvtGenReport( EVTGEN_INFO, "EvtGen" )		EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "TAUOLA negative charge spin propagator PDG id set to "		<< "TAUOLA negative charge spin propagator PDG id set to "
<< _negPropType << endl;		<< m_negPropType << endl;
}		}

if ( _posPropType != 0 ) {		if ( m_posPropType != 0 ) {
EvtGenReport( EVTGEN_INFO, "EvtGen" )		EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "TAUOLA positive charge spin propagator PDG id set to "		<< "TAUOLA positive charge spin propagator PDG id set to "
<< _posPropType << endl;		<< m_posPropType << endl;
}		}

// 3) TauolaHiggsMixingAngle: Specify the mixing angle between the neutral scalar & pseudoscalar Higgs		// 3) TauolaHiggsMixingAngle: Specify the mixing angle between the neutral scalar & pseudoscalar Higgs
// A0/H0; the default mixing angle is pi/4 radians		// A0/H0; the default mixing angle is pi/4 radians
std::string mixString = EvtSymTable::get( "TauolaHiggsMixingAngle", iErr );		std::string mixString = EvtSymTable::get( "TauolaHiggsMixingAngle", iErr );
// If the definition name is not found, get() just returns the first argument string		// If the definition name is not found, get() just returns the first argument string
if ( mixString != "TauolaHiggsMixingAngle" ) {		if ( mixString != "TauolaHiggsMixingAngle" ) {
double mixAngle = std::atof( mixString.c_str() );		double mixAngle = std::atof( mixString.c_str() );
Show All 40 Lines	if ( currentOption != "TauolaCurrentOption" ) {
<< "TAUOLA current option = " << currentOpt << endl;		<< "TAUOLA current option = " << currentOpt << endl;

Tauolapp::Tauola::setNewCurrents( currentOpt );		Tauolapp::Tauola::setNewCurrents( currentOpt );
}		}
}		}

bool EvtTauolaEngine::doDecay( EvtParticle* tauParticle )		bool EvtTauolaEngine::doDecay( EvtParticle* tauParticle )
{		{
if ( _initialised == false ) {		if ( m_initialised == false ) {
this->initialise();		this->initialise();
}		}

if ( tauParticle == 0 ) {		if ( tauParticle == 0 ) {
return false;		return false;
}		}

// Check that we have a tau particle.		// Check that we have a tau particle.
EvtId partId = tauParticle->getId();		EvtId partId = tauParticle->getId();
if ( abs( EvtPDL::getStdHep( partId ) ) != _tauPDG ) {		if ( abs( EvtPDL::getStdHep( partId ) ) != m_tauPDG ) {
return false;		return false;
}		}

int nTauDaug = tauParticle->getNDaug();		int nTauDaug = tauParticle->getNDaug();

// If the number of tau daughters is not zero, then we have already decayed		// If the number of tau daughters is not zero, then we have already decayed
// it using Tauola/another decay algorithm.		// it using Tauola/another decay algorithm.
if ( nTauDaug > 0 ) {		if ( nTauDaug > 0 ) {
▲ Show 20 Lines • Show All 74 Lines • ▼ Show 20 Lines	if ( theParent ) {
if ( theDaughter ) {		if ( theDaughter ) {
GenParticlePtr hepMCDaughter = this->createGenParticle(		GenParticlePtr hepMCDaughter = this->createGenParticle(
theDaughter );		theDaughter );
theVertex->add_particle_out( hepMCDaughter );		theVertex->add_particle_out( hepMCDaughter );

EvtId theId = theDaughter->getId();		EvtId theId = theDaughter->getId();
int PDGInt = EvtPDL::getStdHep( theId );		int PDGInt = EvtPDL::getStdHep( theId );

if ( abs( PDGInt ) == _tauPDG ) {		if ( abs( PDGInt ) == m_tauPDG ) {
// Delete any siblings for the tau particle		// Delete any siblings for the tau particle
if ( theDaughter->getNDaug() > 0 ) {		if ( theDaughter->getNDaug() > 0 ) {
theDaughter->deleteDaughters( false );		theDaughter->deleteDaughters( false );
}		}
tauMap[hepMCDaughter] = theDaughter;		tauMap[hepMCDaughter] = theDaughter;
nTaus++;		nTaus++;
} else {		} else {
// Treat all other particles as "stable"		// Treat all other particles as "stable"
hepMCDaughter->set_status( Tauolapp::TauolaParticle::STABLE );		hepMCDaughter->set_status( Tauolapp::TauolaParticle::STABLE );
}		}

} // theDaughter != 0		} // theDaughter != 0
} // Loop over daughters		} // Loop over daughters

// For the parent particle, artifically set the PDG to a boson with the same 4-momentum		// For the parent particle, artifically set the PDG to a boson with the same 4-momentum
// so that spin correlations are calculated inside Tauola.		// so that spin correlations are calculated inside Tauola.
// This leaves the original parent _EvtParticle_ unchanged		// This leaves the original parent m_EvtParticle_ unchanged
if ( nTaus > 0 && hepMCParent ) {		if ( nTaus > 0 && hepMCParent ) {
int parCharge = EvtPDL::chg3( origParentId ) /		int parCharge = EvtPDL::chg3( origParentId ) /
3; // (3*particle charge)/3 = particle charge		3; // (3*particle charge)/3 = particle charge
if ( parCharge == 0 && _neutPropType != 0 ) {		if ( parCharge == 0 && m_neutPropType != 0 ) {
hepMCParent->set_pdg_id( _neutPropType );		hepMCParent->set_pdg_id( m_neutPropType );
} else if ( parCharge == -1 && _negPropType != 0 ) {		} else if ( parCharge == -1 && m_negPropType != 0 ) {
hepMCParent->set_pdg_id( _negPropType );		hepMCParent->set_pdg_id( m_negPropType );
} else if ( parCharge == 1 && _posPropType != 0 ) {		} else if ( parCharge == 1 && m_posPropType != 0 ) {
hepMCParent->set_pdg_id( _posPropType );		hepMCParent->set_pdg_id( m_posPropType );
}		}
}		}

} else {		} else {
// We only have the one tau particle. Store only this in the map.		// We only have the one tau particle. Store only this in the map.
GenParticlePtr singleTau = this->createGenParticle( tauParticle );		GenParticlePtr singleTau = this->createGenParticle( tauParticle );
theVertex->add_particle_out( singleTau );		theVertex->add_particle_out( singleTau );
tauMap[singleTau] = tauParticle;		tauMap[singleTau] = tauParticle;
Show All 23 Lines
#else		#else
HepMC::GenEvent::particle_iterator eventIter;		HepMC::GenEvent::particle_iterator eventIter;
for ( eventIter = theEvent->particles_begin();		for ( eventIter = theEvent->particles_begin();
eventIter != theEvent->particles_end(); ++eventIter ) {		eventIter != theEvent->particles_end(); ++eventIter ) {
// Check to see if we have a tau particle		// Check to see if we have a tau particle
HepMC::GenParticle* aParticle = ( *eventIter );		HepMC::GenParticle* aParticle = ( *eventIter );
#endif		#endif

if ( aParticle && abs( aParticle->pdg_id() ) == _tauPDG ) {		if ( aParticle && abs( aParticle->pdg_id() ) == m_tauPDG ) {
// Find out what EvtParticle corresponds to the HepMC particle.		// Find out what EvtParticle corresponds to the HepMC particle.
// We need this to create and attach EvtParticle daughters.		// We need this to create and attach EvtParticle daughters.
EvtParticle* tauEvtParticle = tauMap[aParticle];		EvtParticle* tauEvtParticle = tauMap[aParticle];

if ( tauEvtParticle ) {		if ( tauEvtParticle ) {
// Get the tau 4-momentum in the lab (first mother) frame. We need to boost		// Get the tau 4-momentum in the lab (first mother) frame. We need to boost
// all the tau daughters to this frame, such that daug.getP4() is in the tau restframe.		// all the tau daughters to this frame, such that daug.getP4() is in the tau restframe.
EvtVector4R tauP4CM = tauEvtParticle->getP4Lab();		EvtVector4R tauP4CM = tauEvtParticle->getP4Lab();
▲ Show 20 Lines • Show All 101 Lines • Show Last 20 Lines