diff --git a/EvtGenBase/EvtComplex.hh b/EvtGenBase/EvtComplex.hh
index c6e392a..4245f2c 100644
--- a/EvtGenBase/EvtComplex.hh
+++ b/EvtGenBase/EvtComplex.hh
@@ -1,249 +1,250 @@
/***********************************************************************
* 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 . *
***********************************************************************/
#ifndef EVTCOMPLEX_HH
#define EVTCOMPLEX_HH
#include "EvtGenBase/EvtConst.hh"
#include
#include
#include
class EvtComplex {
inline friend EvtComplex operator*( double d, const EvtComplex& c );
inline friend EvtComplex operator*( const EvtComplex& c, double d );
inline friend EvtComplex operator/( const EvtComplex& c, double d );
inline friend EvtComplex operator/( double d, const EvtComplex& c );
inline friend EvtComplex operator*( const EvtComplex& c1,
const EvtComplex& c2 );
inline friend EvtComplex operator/( const EvtComplex& c1,
const EvtComplex& c2 );
inline friend EvtComplex operator+( const EvtComplex& c1,
const EvtComplex& c2 );
inline friend EvtComplex operator-( const EvtComplex& c1,
const EvtComplex& c2 );
inline friend EvtComplex operator-( const EvtComplex& c );
inline friend EvtComplex conj( const EvtComplex& c );
inline friend double abs( const EvtComplex& c );
inline friend double abs2( const EvtComplex& c );
inline friend double arg( const EvtComplex& c );
inline friend double real( const EvtComplex& c );
inline friend double imag( const EvtComplex& c );
inline friend EvtComplex exp( const EvtComplex& c );
friend std::ostream& operator<<( std::ostream& s, const EvtComplex& c );
public:
EvtComplex() : _rpart( 0.0 ), _ipart( 0.0 ) {}
EvtComplex( double rpart, double ipart = 0.0 ) :
_rpart( rpart ), _ipart( ipart )
{
}
EvtComplex( const EvtComplex& c ) : _rpart( c._rpart ), _ipart( c._ipart )
{
}
- EvtComplex( const std::complex& c ) : _rpart( c.real() ), _ipart( c.imag() )
+ EvtComplex( const std::complex& c ) :
+ _rpart( c.real() ), _ipart( c.imag() )
{
- }
+ }
inline EvtComplex& operator*=( double d );
inline EvtComplex& operator/=( double d );
EvtComplex& operator*=( EvtComplex c );
EvtComplex& operator/=( EvtComplex c );
inline EvtComplex& operator=( const EvtComplex& c );
inline EvtComplex& operator+=( const EvtComplex& c );
inline EvtComplex& operator-=( const EvtComplex& c );
inline EvtComplex& operator+=( double d );
inline EvtComplex& operator-=( double d );
inline int operator==( const EvtComplex c );
inline int operator!=( const EvtComplex c );
private:
double _rpart, _ipart;
};
typedef EvtComplex* EvtComplexPtr;
typedef EvtComplexPtr* EvtComplexPtrPtr;
typedef EvtComplexPtrPtr* EvtComplexPtrPtrPtr;
EvtComplex& EvtComplex::operator=( const EvtComplex& c )
{
_rpart = c._rpart;
_ipart = c._ipart;
return *this;
}
EvtComplex& EvtComplex::operator+=( const EvtComplex& c )
{
_rpart += c._rpart;
_ipart += c._ipart;
return *this;
}
EvtComplex& EvtComplex::operator-=( const EvtComplex& c )
{
_rpart -= c._rpart;
_ipart -= c._ipart;
return *this;
}
EvtComplex& EvtComplex::operator+=( double d )
{
_rpart += d;
return *this;
}
EvtComplex& EvtComplex::operator-=( double d )
{
_rpart -= d;
return *this;
}
EvtComplex operator*( double d, const EvtComplex& c )
{
return EvtComplex( c._rpart * d, c._ipart * d );
}
EvtComplex operator*( const EvtComplex& c, double d )
{
return EvtComplex( c._rpart * d, c._ipart * d );
}
EvtComplex operator/( const EvtComplex& c, double d )
{
return EvtComplex( c._rpart / d, c._ipart / d );
}
EvtComplex& EvtComplex::operator*=( double d )
{
_rpart *= d;
_ipart *= d;
return *this;
}
EvtComplex& EvtComplex::operator/=( double d )
{
_rpart /= d;
_ipart /= d;
return *this;
}
EvtComplex operator/( double d, const EvtComplex& c )
{
double Num = d / ( c._rpart * c._rpart + c._ipart * c._ipart );
return EvtComplex( Num * c._rpart, -Num * c._ipart );
}
EvtComplex operator/( const EvtComplex& c1, const EvtComplex& c2 )
{
double inv = 1.0 / ( c2._rpart * c2._rpart + c2._ipart * c2._ipart );
return EvtComplex( inv * ( c1._rpart * c2._rpart + c1._ipart * c2._ipart ),
inv * ( c1._ipart * c2._rpart - c1._rpart * c2._ipart ) );
}
EvtComplex operator*( const EvtComplex& c1, const EvtComplex& c2 )
{
return EvtComplex( c1._rpart * c2._rpart - c1._ipart * c2._ipart,
c1._rpart * c2._ipart + c1._ipart * c2._rpart );
}
EvtComplex operator-( const EvtComplex& c1, const EvtComplex& c2 )
{
return EvtComplex( c1._rpart - c2._rpart, c1._ipart - c2._ipart );
}
EvtComplex operator+( const EvtComplex& c1, const EvtComplex& c2 )
{
return EvtComplex( c1._rpart + c2._rpart, c1._ipart + c2._ipart );
}
int EvtComplex::operator==( const EvtComplex c )
{
return _rpart == c._rpart && _ipart == c._ipart;
}
int EvtComplex::operator!=( const EvtComplex c )
{
return _rpart != c._rpart || _ipart != c._ipart;
}
EvtComplex operator-( const EvtComplex& c )
{
return EvtComplex( -c._rpart, -c._ipart );
}
EvtComplex conj( const EvtComplex& c )
{
return EvtComplex( c._rpart, -c._ipart );
}
double abs( const EvtComplex& c )
{
double c2 = c._rpart * c._rpart + c._ipart * c._ipart;
if ( c2 <= 0.0 )
return 0.0;
return sqrt( c2 );
}
double abs2( const EvtComplex& c )
{
return c._rpart * c._rpart + c._ipart * c._ipart;
}
double arg( const EvtComplex& c )
{
if ( ( c._rpart == 0 ) && ( c._ipart == 0 ) ) {
return 0.0;
}
if ( c._rpart == 0 ) {
if ( c._ipart > 0 ) {
return EvtConst::pi / 2;
} else {
return -EvtConst::pi / 2;
}
} else {
return atan2( c._ipart, c._rpart );
}
}
double real( const EvtComplex& c )
{
return c._rpart;
}
double imag( const EvtComplex& c )
{
return c._ipart;
}
EvtComplex exp( const EvtComplex& c )
{
return exp( c._rpart ) * EvtComplex( cos( c._ipart ), sin( c._ipart ) );
}
#endif
diff --git a/src/EvtGenExternal/EvtExternalGenList.cpp b/src/EvtGenExternal/EvtExternalGenList.cpp
index 003706d..5414bb7 100644
--- a/src/EvtGenExternal/EvtExternalGenList.cpp
+++ b/src/EvtGenExternal/EvtExternalGenList.cpp
@@ -1,78 +1,82 @@
/***********************************************************************
* 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 . *
***********************************************************************/
#include "EvtGenExternal/EvtExternalGenList.hh"
#include "EvtGenExternal/EvtExternalGenFactory.hh"
#include "EvtGenExternal/EvtPHOTOS.hh"
#include "EvtGenExternal/EvtPythia.hh"
#include "EvtGenExternal/EvtTauola.hh"
+#include "EvtGenExternal/EvtHmePythia.hh"
EvtExternalGenList::EvtExternalGenList( bool convertPythiaCodes,
std::string pythiaXmlDir,
std::string photonType,
bool 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
// xmldoc Pythia directory
char* pythiaDataDir = getenv( "PYTHIA8DATA" );
if ( pythiaDataDir != 0 ) {
pythiaXmlDir = pythiaDataDir;
}
}
extFactory->definePythiaGenerator( pythiaXmlDir, convertPythiaCodes,
useEvtGenRandom );
extFactory->defineTauolaGenerator( useEvtGenRandom );
}
EvtExternalGenList::~EvtExternalGenList()
{
}
EvtAbsRadCorr* EvtExternalGenList::getPhotosModel()
{
// Define the Photos model, which uses the EvtPhotosEngine class.
EvtPHOTOS* photosModel = new EvtPHOTOS();
return photosModel;
}
std::list EvtExternalGenList::getListOfModels()
{
// Create the Pythia and Tauola models, which use their own engine classes.
EvtPythia* pythiaModel = new EvtPythia();
EvtTauola* tauolaModel = new EvtTauola();
+ // Pythia helicity matrix element decay model, e.g. for tau decays
+ EvtHmePythia* hmePythiaModel = new EvtHmePythia();
std::list extraModels;
extraModels.push_back( pythiaModel );
extraModels.push_back( tauolaModel );
+ extraModels.push_back( hmePythiaModel );
// Return the list of models
return extraModels;
}
diff --git a/src/EvtGenExternal/EvtHmePythia.cpp b/src/EvtGenExternal/EvtHmePythia.cpp
index 2148415..55a9013 100644
--- a/src/EvtGenExternal/EvtHmePythia.cpp
+++ b/src/EvtGenExternal/EvtHmePythia.cpp
@@ -1,77 +1,90 @@
/***********************************************************************
* 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 . *
***********************************************************************/
#include "EvtGenExternal/EvtHmePythia.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenExternal/EvtExternalGenFactory.hh"
#include "EvtGenExternal/EvtPythiaEngine.hh"
// Name
std::string EvtHmePythia::getName()
{
return "HMEPYTHIA";
}
// Clone the model
EvtDecayBase* EvtHmePythia::clone()
{
return new EvtHmePythia();
}
// Initialize the model
void EvtHmePythia::init()
{
// Retrieve the Pythia engine
if ( !_pythiaEngine ) {
_pythiaEngine = dynamic_cast(
EvtExternalGenFactory::getInstance()->getGenerator(
EvtExternalGenFactory::PythiaGenId ) );
}
- // Check the matrix element is specified and is available
+ // Check the matrix element is specified and is available.
+ // Allowed values are (all implicitly include the tau neutrino):
+ // 1521 = HMETau2Meson
+ // 1531 = HMETau2TwoLeptons
+ // 1532 = HMETau2TwoMesonsViaVector
+ // 1533 = HMETau2TwoMesonsViaVectorScalar
+ // 1541 = HMETau2ThreePions
+ // 1542 = HMETau2ThreeMesonsWithKaons
+ // 1543 = HMETau2ThreeMesonsGeneric
+ // 1544 = HMETau2TwoPionsGamma
+ // 1551 = HMETau2FourPions
+ // 1561 = HMETau2FivePions
+
checkNArg( 1 );
_hmeMode = getArg( 0 );
if ( _pythiaEngine->hmeAmplitudes( _hmeMode ) == 0 ) {
EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< "Helicity matrix element mode " << _hmeMode
<< " is not available." << std::endl;
}
// Check the number of daughters matches the matrix element
+ // e.g. tau- -> e- anti-nu_e nu_tau, HMETau2TwoLeptons (1531), 3 daughters
checkNDaug( ( _hmeMode % 100 ) / 10 );
}
void EvtHmePythia::initProbMax()
{
setProbMax( _pythiaEngine->hmeAmplitudes( _hmeMode ) );
}
void EvtHmePythia::decay( EvtParticle* p )
{
// Initialize the phase-space
p->initializePhaseSpace( getNDaug(), getDaugs() );
// Call the Pythia helicity matrix element
if ( _pythiaEngine ) {
_pythiaEngine->hmeAmplitudes( _hmeMode, p, this );
}
}