diff --git a/src/EvtGenBase/EvtGenKine.cpp b/src/EvtGenBase/EvtGenKine.cpp index c80455c..95c94f8 100644 --- a/src/EvtGenBase/EvtGenKine.cpp +++ b/src/EvtGenBase/EvtGenKine.cpp @@ -1,400 +1,401 @@ /*********************************************************************** * 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 "EvtGenBase/EvtGenKine.hh" #include "EvtGenBase/EvtConst.hh" +#include "EvtGenBase/EvtParticle.hh" #include "EvtGenBase/EvtPatches.hh" #include "EvtGenBase/EvtRandom.hh" #include "EvtGenBase/EvtReport.hh" #include "EvtGenBase/EvtVector4R.hh" -#include "EvtGenBase/EvtParticle.hh" +#include #include -#include + using std::endl; double EvtPawt( double a, double b, double c ) { double temp = ( a * a - ( b + c ) * ( b + c ) ) * ( a * a - ( b - c ) * ( b - c ) ); if ( temp <= 0 ) { return 0.0; } return sqrt( temp ) / ( 2.0 * a ); } double EvtGenKine::PhaseSpace( int ndaug, double mass[30], EvtVector4R p4[30], double mp ) // N body phase space routine. Send parent with // daughters already defined ( Number and masses ) // Returns four vectors in parent frame. { double energy, p3, alpha, beta; if ( ndaug == 1 ) { p4[0].set( mass[0], 0.0, 0.0, 0.0 ); return 1.0; } if ( ndaug == 2 ) { //Two body phase space energy = ( mp * mp + mass[0] * mass[0] - mass[1] * mass[1] ) / ( 2.0 * mp ); p3 = 0.0; if ( energy > mass[0] ) { p3 = sqrt( energy * energy - mass[0] * mass[0] ); } p4[0].set( energy, 0.0, 0.0, p3 ); energy = mp - energy; p3 = -1.0 * p3; p4[1].set( energy, 0.0, 0.0, p3 ); //Now rotate four vectors. alpha = EvtRandom::Flat( EvtConst::twoPi ); beta = acos( EvtRandom::Flat( -1.0, 1.0 ) ); p4[0].applyRotateEuler( alpha, beta, -alpha ); p4[1].applyRotateEuler( alpha, beta, -alpha ); return 1.0; } if ( ndaug != 2 ) { double wtmax = 0.0; double pm[5][30], pmin, pmax, psum, rnd[30]; double ran, wt, pa, costh, sinth, phi, p[4][30], be[4], bep, temp; int i, il, ilr, i1, il1u, il1, il2r, ilu; int il2 = 0; for ( i = 0; i < ndaug; i++ ) { pm[4][i] = 0.0; rnd[i] = 0.0; } pm[0][0] = mp; pm[1][0] = 0.0; pm[2][0] = 0.0; pm[3][0] = 0.0; pm[4][0] = mp; psum = 0.0; for ( i = 1; i < ndaug + 1; i++ ) { psum = psum + mass[i - 1]; } pm[4][ndaug - 1] = mass[ndaug - 1]; switch ( ndaug ) { case 1: wtmax = 1.0 / 16.0; break; case 2: wtmax = 1.0 / 150.0; break; case 3: wtmax = 1.0 / 2.0; break; case 4: wtmax = 1.0 / 5.0; break; case 5: wtmax = 1.0 / 15.0; break; case 6: wtmax = 1.0 / 15.0; break; case 7: wtmax = 1.0 / 15.0; break; case 8: wtmax = 1.0 / 15.0; break; case 9: wtmax = 1.0 / 15.0; break; case 10: wtmax = 1.0 / 15.0; break; case 11: wtmax = 1.0 / 15.0; break; case 12: wtmax = 1.0 / 15.0; break; case 13: wtmax = 1.0 / 15.0; break; case 14: wtmax = 1.0 / 15.0; break; case 15: wtmax = 1.0 / 15.0; break; default: EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << "too many daughters for phase space..." << ndaug << " " << mp << endl; ; break; } pmax = mp - psum + mass[ndaug - 1]; pmin = 0.0; for ( ilr = 2; ilr < ndaug + 1; ilr++ ) { il = ndaug + 1 - ilr; pmax = pmax + mass[il - 1]; pmin = pmin + mass[il + 1 - 1]; wtmax = wtmax * EvtPawt( pmax, pmin, mass[il - 1] ); } do { rnd[0] = 1.0; il1u = ndaug - 1; for ( il1 = 2; il1 < il1u + 1; il1++ ) { ran = EvtRandom::Flat(); for ( il2r = 2; il2r < il1 + 1; il2r++ ) { il2 = il1 + 1 - il2r; if ( ran <= rnd[il2 - 1] ) goto two39; rnd[il2 + 1 - 1] = rnd[il2 - 1]; } two39: rnd[il2 + 1 - 1] = ran; } rnd[ndaug - 1] = 0.0; wt = 1.0; for ( ilr = 2; ilr < ndaug + 1; ilr++ ) { il = ndaug + 1 - ilr; pm[4][il - 1] = pm[4][il + 1 - 1] + mass[il - 1] + ( rnd[il - 1] - rnd[il + 1 - 1] ) * ( mp - psum ); wt = wt * EvtPawt( pm[4][il - 1], pm[4][il + 1 - 1], mass[il - 1] ); } if ( wt > wtmax ) { EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << "wtmax to small in EvtPhaseSpace with " << ndaug << " daughters" << endl; ; } } while ( wt < EvtRandom::Flat( wtmax ) ); ilu = ndaug - 1; for ( il = 1; il < ilu + 1; il++ ) { pa = EvtPawt( pm[4][il - 1], pm[4][il + 1 - 1], mass[il - 1] ); costh = EvtRandom::Flat( -1.0, 1.0 ); sinth = sqrt( 1.0 - costh * costh ); phi = EvtRandom::Flat( EvtConst::twoPi ); p[1][il - 1] = pa * sinth * cos( phi ); p[2][il - 1] = pa * sinth * sin( phi ); p[3][il - 1] = pa * costh; pm[1][il + 1 - 1] = -p[1][il - 1]; pm[2][il + 1 - 1] = -p[2][il - 1]; pm[3][il + 1 - 1] = -p[3][il - 1]; p[0][il - 1] = sqrt( pa * pa + mass[il - 1] * mass[il - 1] ); pm[0][il + 1 - 1] = sqrt( pa * pa + pm[4][il + 1 - 1] * pm[4][il + 1 - 1] ); } p[0][ndaug - 1] = pm[0][ndaug - 1]; p[1][ndaug - 1] = pm[1][ndaug - 1]; p[2][ndaug - 1] = pm[2][ndaug - 1]; p[3][ndaug - 1] = pm[3][ndaug - 1]; for ( ilr = 2; ilr < ndaug + 1; ilr++ ) { il = ndaug + 1 - ilr; be[0] = pm[0][il - 1] / pm[4][il - 1]; be[1] = pm[1][il - 1] / pm[4][il - 1]; be[2] = pm[2][il - 1] / pm[4][il - 1]; be[3] = pm[3][il - 1] / pm[4][il - 1]; for ( i1 = il; i1 < ndaug + 1; i1++ ) { bep = be[1] * p[1][i1 - 1] + be[2] * p[2][i1 - 1] + be[3] * p[3][i1 - 1] + be[0] * p[0][i1 - 1]; temp = ( p[0][i1 - 1] + bep ) / ( be[0] + 1.0 ); p[1][i1 - 1] = p[1][i1 - 1] + temp * be[1]; p[2][i1 - 1] = p[2][i1 - 1] + temp * be[2]; p[3][i1 - 1] = p[3][i1 - 1] + temp * be[3]; p[0][i1 - 1] = bep; } } for ( ilr = 0; ilr < ndaug; ilr++ ) { p4[ilr].set( p[0][ilr], p[1][ilr], p[2][ilr], p[3][ilr] ); } return 1.0; } return 1.0; } double EvtGenKine::PhaseSpacePole( double M, double m1, double m2, double m3, double a, EvtVector4R p4[10] ) // generate kinematics for 3 body decays, pole for the m1,m2 mass. { //f1 = 1 (phasespace) //f2 = a*(1/m12sq)^2 double m12sqmax = ( M - m3 ) * ( M - m3 ); double m12sqmin = ( m1 + m2 ) * ( m1 + m2 ); double m13sqmax = ( M - m2 ) * ( M - m2 ); double m13sqmin = ( m1 + m3 ) * ( m1 + m3 ); double v1 = ( m12sqmax - m12sqmin ) * ( m13sqmax - m13sqmin ); double v2 = a * ( 1.0 / m12sqmin - 1.0 / m12sqmax ) * ( m13sqmax - m13sqmin ); double m12sq, m13sq; double r = v1 / ( v1 + v2 ); double m13min, m13max; do { m13sq = EvtRandom::Flat( m13sqmin, m13sqmax ); if ( r > EvtRandom::Flat() ) { m12sq = EvtRandom::Flat( m12sqmin, m12sqmax ); } else { m12sq = 1.0 / ( 1.0 / m12sqmin - EvtRandom::Flat() * ( 1.0 / m12sqmin - 1.0 / m12sqmax ) ); } //kinematically allowed? double E3star = ( M * M - m12sq - m3 * m3 ) / sqrt( 4 * m12sq ); double E1star = ( m12sq + m1 * m1 - m2 * m2 ) / sqrt( 4 * m12sq ); double p3star = sqrt( E3star * E3star - m3 * m3 ); double p1star = sqrt( E1star * E1star - m1 * m1 ); m13max = ( E3star + E1star ) * ( E3star + E1star ) - ( p3star - p1star ) * ( p3star - p1star ); m13min = ( E3star + E1star ) * ( E3star + E1star ) - ( p3star + p1star ) * ( p3star + p1star ); } while ( m13sq < m13min || m13sq > m13max ); double E2 = ( M * M + m2 * m2 - m13sq ) / ( 2.0 * M ); double E3 = ( M * M + m3 * m3 - m12sq ) / ( 2.0 * M ); double E1 = M - E2 - E3; double p1mom = sqrt( E1 * E1 - m1 * m1 ); double p3mom = sqrt( E3 * E3 - m3 * m3 ); double cost13 = ( 2.0 * E1 * E3 + m1 * m1 + m3 * m3 - m13sq ) / ( 2.0 * p1mom * p3mom ); //EvtGenReport(EVTGEN_INFO,"EvtGen") << m13sq << endl; //EvtGenReport(EVTGEN_INFO,"EvtGen") << m12sq << endl; //EvtGenReport(EVTGEN_INFO,"EvtGen") << E1 << endl; //EvtGenReport(EVTGEN_INFO,"EvtGen") << E2 << endl; //EvtGenReport(EVTGEN_INFO,"EvtGen") << E3 << endl; //EvtGenReport(EVTGEN_INFO,"EvtGen") << p1mom << endl; //EvtGenReport(EVTGEN_INFO,"EvtGen") << p3mom << endl; //EvtGenReport(EVTGEN_INFO,"EvtGen") << cost13 << endl; p4[2].set( E3, 0.0, 0.0, p3mom ); p4[0].set( E1, p1mom * sqrt( 1.0 - cost13 * cost13 ), 0.0, p1mom * cost13 ); p4[1].set( E2, -p1mom * sqrt( 1.0 - cost13 * cost13 ), 0.0, -p1mom * cost13 - p3mom ); //EvtGenReport(EVTGEN_INFO,"EvtGen") << "p4:"<mass(); EvtParticle* daug1 = p->getDaug( 0 ); EvtParticle* daug2 = p->getDaug( 1 ); EvtParticle* daug3 = p->getDaug( 2 ); const double mDaug1 = daug1->mass(); const double mDaug2 = daug2->mass(); const double mDaug3 = daug3->mass(); const double mParentSq{ mParent * mParent }; const double mDaug1Sq{ mDaug1 * mDaug1 }; const double mDaug2Sq{ mDaug2 * mDaug2 }; const double mDaug3Sq{ mDaug3 * mDaug3 }; const double invMParent{ 1. / mParent }; const double En1 = 0.5 * ( mParentSq + mDaug1Sq - m23Sq ) * invMParent; const double En3 = 0.5 * ( mParentSq + mDaug3Sq - m12Sq ) * invMParent; const double En2 = mParent - En1 - En3; const double p1mag = std::sqrt( En1 * En1 - mDaug1Sq ); const double p2mag = std::sqrt( En2 * En2 - mDaug2Sq ); double cosPhi = 0.5 * ( mDaug1Sq + mDaug2Sq + 2 * En1 * En2 - m12Sq ) / ( p1mag * p2mag ); double sinPhi = std::sqrt( 1 - cosPhi * cosPhi ); if ( EvtRandom::Flat( 0., 1. ) > 0.5 ) { sinPhi *= -1; } const double p2x = p2mag * cosPhi; const double p2y = p2mag * sinPhi; const double p3x = -p1mag - p2x; const double p3y = -p2y; // Construct 4-momenta and rotate them randomly in space EvtVector4R p1( En1, p1mag, 0., 0. ); EvtVector4R p2( En2, p2x, p2y, 0. ); EvtVector4R p3( En3, p3x, p3y, 0. ); const double euler1 = EvtRandom::Flat( 0., EvtConst::twoPi ); const double euler2 = std::acos( EvtRandom::Flat( -1.0, 1.0 ) ); const double euler3 = EvtRandom::Flat( 0., EvtConst::twoPi ); p1.applyRotateEuler( euler1, euler2, euler3 ); p2.applyRotateEuler( euler1, euler2, euler3 ); p3.applyRotateEuler( euler1, euler2, euler3 ); // set momenta for daughters daug1->init( daug1->getId(), p1 ); daug2->init( daug2->getId(), p2 ); daug3->init( daug3->getId(), p3 ); return; } diff --git a/src/EvtGenExternal/EvtExternalGenFactory.cpp b/src/EvtGenExternal/EvtExternalGenFactory.cpp index ce9e640..b96fdae 100644 --- a/src/EvtGenExternal/EvtExternalGenFactory.cpp +++ b/src/EvtGenExternal/EvtExternalGenFactory.cpp @@ -1,159 +1,165 @@ /*********************************************************************** * 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/EvtExternalGenFactory.hh" #include "EvtGenBase/EvtPatches.hh" #include "EvtGenBase/EvtReport.hh" #ifdef EVTGEN_PYTHIA #include "EvtGenExternal/EvtPythiaEngine.hh" #endif #ifdef EVTGEN_PHOTOS #include "EvtGenExternal/EvtPhotosEngine.hh" #endif #ifdef EVTGEN_TAUOLA #include "EvtGenExternal/EvtTauolaEngine.hh" #endif #include using std::endl; EvtExternalGenFactory::EvtExternalGenFactory() { _extGenMap.clear(); } EvtExternalGenFactory::~EvtExternalGenFactory() { ExtGenMap::iterator iter; for ( iter = _extGenMap.begin(); iter != _extGenMap.end(); ++iter ) { EvtAbsExternalGen* theGenerator = iter->second; delete theGenerator; } _extGenMap.clear(); } EvtExternalGenFactory* EvtExternalGenFactory::getInstance() { static EvtExternalGenFactory* theFactory = 0; if ( theFactory == 0 ) { theFactory = new EvtExternalGenFactory(); } return theFactory; } +// Only define the generator if we have the external ifdef variable set +#ifdef EVTGEN_PYTHIA void EvtExternalGenFactory::definePythiaGenerator( std::string xmlDir, bool convertPhysCodes, bool useEvtGenRandom ) { - // Only define the generator if we have the external ifdef variable set -#ifdef EVTGEN_PYTHIA - int genId = EvtExternalGenFactory::PythiaGenId; EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Defining EvtPythiaEngine: data tables defined in " << xmlDir << endl; if ( convertPhysCodes == true ) { EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Pythia 6 codes in decay files will be converted to Pythia 8 codes" << endl; } else { EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Pythia 8 codes need to be used in decay files" << endl; } if ( useEvtGenRandom == true ) { EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Using EvtGen random engine for Pythia 8 as well" << endl; } EvtAbsExternalGen* pythiaGenerator = new EvtPythiaEngine( xmlDir, convertPhysCodes, useEvtGenRandom ); _extGenMap[genId] = pythiaGenerator; - -#endif } +#else +void EvtExternalGenFactory::definePythiaGenerator( std::string, bool, bool ) +{ +} +#endif +#ifdef EVTGEN_PHOTOS void EvtExternalGenFactory::definePhotosGenerator( std::string photonType, bool useEvtGenRandom ) { -#ifdef EVTGEN_PHOTOS - int genId = EvtExternalGenFactory::PhotosGenId; EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Defining EvtPhotosEngine using photonType = " << photonType << endl; EvtAbsExternalGen* photosGenerator = new EvtPhotosEngine( photonType, useEvtGenRandom ); _extGenMap[genId] = photosGenerator; - -#endif } +#else +void EvtExternalGenFactory::definePhotosGenerator( std::string, bool ) +{ +} +#endif +#ifdef EVTGEN_TAUOLA void EvtExternalGenFactory::defineTauolaGenerator( bool useEvtGenRandom ) { -#ifdef EVTGEN_TAUOLA - int genId = EvtExternalGenFactory::TauolaGenId; EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Defining EvtTauolaEngine." << endl; EvtAbsExternalGen* tauolaGenerator = new EvtTauolaEngine( useEvtGenRandom ); _extGenMap[genId] = tauolaGenerator; - -#endif } +#else +void EvtExternalGenFactory::defineTauolaGenerator( bool ) +{ +} +#endif EvtAbsExternalGen* EvtExternalGenFactory::getGenerator( int genId ) { EvtAbsExternalGen* theGenerator( 0 ); ExtGenMap::iterator iter; if ( ( iter = _extGenMap.find( genId ) ) != _extGenMap.end() ) { // Retrieve the external generator engine theGenerator = iter->second; } else { EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "EvtAbsExternalGen::getGenerator: could not find generator for genId = " << genId << endl; } return theGenerator; } void EvtExternalGenFactory::initialiseAllGenerators() { ExtGenMap::iterator iter; for ( iter = _extGenMap.begin(); iter != _extGenMap.end(); ++iter ) { EvtAbsExternalGen* theGenerator = iter->second; if ( theGenerator != 0 ) { theGenerator->initialise(); } } } diff --git a/src/EvtGenModels/EvtFlatSqDalitz.cpp b/src/EvtGenModels/EvtFlatSqDalitz.cpp index be30e3d..cb46356 100644 --- a/src/EvtGenModels/EvtFlatSqDalitz.cpp +++ b/src/EvtGenModels/EvtFlatSqDalitz.cpp @@ -1,115 +1,112 @@ /*********************************************************************** * 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 "EvtGenModels/EvtFlatSqDalitz.hh" #include "EvtGenBase/EvtConst.hh" #include "EvtGenBase/EvtGenKine.hh" #include "EvtGenBase/EvtPDL.hh" #include "EvtGenBase/EvtParticle.hh" #include "EvtGenBase/EvtPatches.hh" #include "EvtGenBase/EvtRandom.hh" #include "EvtGenBase/EvtReport.hh" -#include -#include -#include +#include #include -using std::fstream; EvtFlatSqDalitz::~EvtFlatSqDalitz() { } std::string EvtFlatSqDalitz::getName() { return "FLATSQDALITZ"; } EvtDecayBase* EvtFlatSqDalitz::clone() { return new EvtFlatSqDalitz; } void EvtFlatSqDalitz::initProbMax() { noProbMax(); } void EvtFlatSqDalitz::init() { //check there are 3 daughters checkNDaug( 3 ); // check that there are 0 arguments checkNArg( 0, 2, 4 ); if ( getNArg() > 0 ) { m_mPrimeMin = getArg( 0 ); m_mPrimeMax = getArg( 1 ); } if ( getNArg() > 2 ) { m_thetaPrimeMin = getArg( 2 ); m_thetaPrimeMax = getArg( 3 ); } } void EvtFlatSqDalitz::decay( EvtParticle* p ) { p->makeDaughters( getNDaug(), getDaugs() ); p->generateMassTree(); const double mParent = p->mass(); EvtParticle* daug1 = p->getDaug( 0 ); EvtParticle* daug2 = p->getDaug( 1 ); EvtParticle* daug3 = p->getDaug( 2 ); const double mDaug1 = daug1->mass(); const double mDaug2 = daug2->mass(); const double mDaug3 = daug3->mass(); const double mParentSq = mParent * mParent; const double mDaug1Sq = mDaug1 * mDaug1; const double mDaug2Sq = mDaug2 * mDaug2; const double mDaug3Sq = mDaug3 * mDaug3; // Generate m' and theta' const double mPrime = EvtRandom::Flat( m_mPrimeMin, m_mPrimeMax ); const double thetaPrime = EvtRandom::Flat( m_thetaPrimeMin, m_thetaPrimeMax ); // calculate m12 and m23 const double m12 = 0.5 * ( std::cos( mPrime * EvtConst::pi ) + 1 ) * ( mParent - ( mDaug1 + mDaug2 + mDaug3 ) ) + mDaug1 + mDaug2; const double m12Sq = m12 * m12; const double en1 = ( m12Sq - mDaug2Sq + mDaug1Sq ) / ( 2. * m12 ); const double en3 = ( mParentSq - m12Sq - mDaug3Sq ) / ( 2. * m12 ); const double p1 = std::sqrt( en1 * en1 - mDaug1Sq ); const double p3 = std::sqrt( en3 * en3 - mDaug3Sq ); const double m13Sq = mDaug1Sq + mDaug3Sq + 2.0 * ( en1 * en3 - p1 * p3 * std::cos( EvtConst::pi * thetaPrime ) ); const double m23Sq = mParentSq - m12Sq - m13Sq + mDaug1Sq + mDaug2Sq + mDaug3Sq; // Turn m12 and m23 into momenta EvtGenKine::ThreeBodyKine( m12Sq, m23Sq, p ); return; }