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 ); } }