Skip to content

Introduction of massive RAMBO in standalone_mg7. #5

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
stloufra merged 5 commits into from
Jun 15, 2026
Merged

Introduction of massive RAMBO in standalone_mg7. #5

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
0664e6c
Umami addition of get meta masses
stloufra Jun 11, 2026
b462399
Addition of RAMBO respecting mass
stloufra Jun 11, 2026
52a6d94
RamboSamplingKernels classic_rambo implementation
stloufra Jun 11, 2026
a111d96
Massive RAMBO implementation in check_sa
stloufra Jun 11, 2026
6f599e8
PR suggestion implementation
stloufra Jun 15, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
120 changes: 91 additions & 29 deletions madgraph/iolibs/template_files/madmatrix/RamboSamplingKernels.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,8 @@
#include "MemoryAccessRandomNumbers.h"
#include "MemoryAccessWeights.h"
#include "MemoryBuffers.h"
#include "rambo.h" // inline implementation of RAMBO algorithms and kernels
#include "rambo.h" // inline classic (massive) RAMBO, ported from standalone_cpp
#include "massless_rambo.h" // inline implementation of massless RAMBO algorithms and kernels

#include <sstream>

Expand All @@ -23,27 +24,27 @@ namespace mg5amcCpu
{
//--------------------------------------------------------------------------

RamboSamplingKernelHost::RamboSamplingKernelHost( const fptype energy, // input: energy
MasslessRamboSamplingKernelHost::MasslessRamboSamplingKernelHost( const fptype energy, // input: energy
const BufferRndNumMomenta& rndmom, // input: random numbers in [0,1]
BufferMomenta& momenta, // output: momenta
BufferWeights& weights, // output: weights
const size_t nevt )
: SamplingKernelBase( energy, rndmom, momenta, weights )
, NumberOfEvents( nevt )
{
if( m_rndmom.isOnDevice() ) throw std::runtime_error( "RamboSamplingKernelHost: rndmom must be a host array" );
if( m_momenta.isOnDevice() ) throw std::runtime_error( "RamboSamplingKernelHost: momenta must be a host array" );
if( m_weights.isOnDevice() ) throw std::runtime_error( "RamboSamplingKernelHost: weights must be a host array" );
if( this->nevt() != m_rndmom.nevt() ) throw std::runtime_error( "RamboSamplingKernelHost: nevt mismatch with rndmom" );
if( this->nevt() != m_momenta.nevt() ) throw std::runtime_error( "RamboSamplingKernelHost: nevt mismatch with momenta" );
if( this->nevt() != m_weights.nevt() ) throw std::runtime_error( "RamboSamplingKernelHost: nevt mismatch with weights" );
if( m_rndmom.isOnDevice() ) throw std::runtime_error( "MasslessRamboSamplingKernelHost: rndmom must be a host array" );
if( m_momenta.isOnDevice() ) throw std::runtime_error( "MasslessRamboSamplingKernelHost: momenta must be a host array" );
if( m_weights.isOnDevice() ) throw std::runtime_error( "MasslessRamboSamplingKernelHost: weights must be a host array" );
if( this->nevt() != m_rndmom.nevt() ) throw std::runtime_error( "MasslessRamboSamplingKernelHost: nevt mismatch with rndmom" );
if( this->nevt() != m_momenta.nevt() ) throw std::runtime_error( "MasslessRamboSamplingKernelHost: nevt mismatch with momenta" );
if( this->nevt() != m_weights.nevt() ) throw std::runtime_error( "MasslessRamboSamplingKernelHost: nevt mismatch with weights" );
// Sanity checks for memory access (momenta buffer)
constexpr int neppM = MemoryAccessMomenta::neppM; // AOSOA layout
static_assert( ispoweroftwo( neppM ), "neppM is not a power of 2" );
if( nevt % neppM != 0 )
{
std::ostringstream sstr;
sstr << "RamboSamplingKernelHost: nevt should be a multiple of neppM=" << neppM;
sstr << "MasslessRamboSamplingKernelHost: nevt should be a multiple of neppM=" << neppM;
throw std::runtime_error( sstr.str() );
}
// Sanity checks for memory access (random number buffer)
Expand All @@ -52,17 +53,17 @@ namespace mg5amcCpu
if( nevt % neppR != 0 )
{
std::ostringstream sstr;
sstr << "RamboSamplingKernelHost: nevt should be a multiple of neppR=" << neppR;
sstr << "MasslessRamboSamplingKernelHost: nevt should be a multiple of neppR=" << neppR;
throw std::runtime_error( sstr.str() );
}
}

//--------------------------------------------------------------------------

void
RamboSamplingKernelHost::getMomentaInitial()
MasslessRamboSamplingKernelHost::getMomentaInitial()
{
constexpr auto getMomentaInitial = ramboGetMomentaInitial<HostAccessMomenta>;
constexpr auto getMomentaInitial = massless_rambo::ramboGetMomentaInitial<HostAccessMomenta>;
// ** START LOOP ON IEVT **
for( size_t ievt = 0; ievt < nevt(); ++ievt )
{
Expand All @@ -76,9 +77,9 @@ namespace mg5amcCpu
//--------------------------------------------------------------------------

void
RamboSamplingKernelHost::getMomentaFinal()
MasslessRamboSamplingKernelHost::getMomentaFinal()
{
constexpr auto getMomentaFinal = ramboGetMomentaFinal<HostAccessRandomNumbers, HostAccessMomenta, HostAccessWeights>;
constexpr auto getMomentaFinal = massless_rambo::ramboGetMomentaFinal<HostAccessRandomNumbers, HostAccessMomenta, HostAccessWeights>;
// ** START LOOP ON IEVT **
for( size_t ievt = 0; ievt < nevt(); ++ievt )
{
Expand All @@ -93,8 +94,69 @@ namespace mg5amcCpu

//--------------------------------------------------------------------------

RamboSamplingKernelHost::RamboSamplingKernelHost( const fptype energy, // input: energy
const BufferRndNumMomenta& rndmom, // input: random [0,1] UNUSED
const std::vector<fptype>& masses, // input: external-leg masses
const int ninitial, // input: #initial-state particles
const size_t nevt, // input: #events
BufferMomenta& momenta, // output: momenta
BufferWeights& weights ) // output: weights
: SamplingKernelBase( energy, rndmom, momenta, weights )
, NumberOfEvents( nevt )
, m_masses( masses.begin(), masses.end() )
, m_ninitial( ninitial )
{
if( m_momenta.isOnDevice() ) throw std::runtime_error( "RamboSamplingKernelHost: momenta must be a host array" );
if( m_weights.isOnDevice() ) throw std::runtime_error( "RamboSamplingKernelHost: weights must be a host array" );
if( this->nevt() != m_momenta.nevt() ) throw std::runtime_error( "RamboSamplingKernelHost: nevt mismatch with momenta" );
if( this->nevt() != m_weights.nevt() ) throw std::runtime_error( "RamboSamplingKernelHost: nevt mismatch with weights" );

constexpr int neppM = MemoryAccessMomenta::neppM; // AOSOA layout

static_assert( ispoweroftwo( neppM ), "neppM is not a power of 2" );
if( nevt % neppM != 0 )
{
std::ostringstream sstr;
sstr << "RamboSamplingKernelHost: nevt should be a multiple of neppM=" << neppM;
throw std::runtime_error( sstr.str() );
}
}

//--------------------------------------------------------------------------

void
RamboSamplingKernelHost::getMomentaInitial()
{
// NOOP
}

//--------------------------------------------------------------------------

void
RamboSamplingKernelHost::getMomentaFinal()
{
const int npar = (int)m_masses.size();
// ** START LOOP ON IEVT **
for( size_t ievt = 0; ievt < nevt(); ++ievt )
{
// Clas. RAMBO returns [E,px,py,pz] vector per ex. particle
// own RNG, intial final once
// For reproducibility betwn fptype = FP32/FP64 generation in FP64
double wgt = 0.;
const std::vector<std::vector<double>> point =
rambo::get_momenta( m_ninitial, (double)m_energy, m_masses, wgt );
for( int ipar = 0; ipar < npar; ++ipar )
for( int ip4 = 0; ip4 < 4; ++ip4 )
MemoryAccessMomenta::ieventAccessIp4Ipar( m_momenta.data(), ievt, ip4, ipar ) = (fptype)point[ipar][ip4];
MemoryAccessWeights::ieventAccess( m_weights.data(), ievt ) = (fptype)wgt;
}
// ** END LOOP ON IEVT **
}

//--------------------------------------------------------------------------

#ifdef MGONGPUCPP_GPUIMPL
RamboSamplingKernelDevice::RamboSamplingKernelDevice( const fptype energy, // input: energy
MasslessRamboSamplingKernelDevice::MasslessRamboSamplingKernelDevice( const fptype energy, // input: energy
const BufferRndNumMomenta& rndmom, // input: random numbers in [0,1]
BufferMomenta& momenta, // output: momenta
BufferWeights& weights, // output: weights
Expand All @@ -105,21 +167,21 @@ namespace mg5amcCpu
, m_gpublocks( gpublocks )
, m_gputhreads( gputhreads )
{
if( !m_rndmom.isOnDevice() ) throw std::runtime_error( "RamboSamplingKernelDevice: rndmom must be a device array" );
if( !m_momenta.isOnDevice() ) throw std::runtime_error( "RamboSamplingKernelDevice: momenta must be a device array" );
if( !m_weights.isOnDevice() ) throw std::runtime_error( "RamboSamplingKernelDevice: weights must be a device array" );
if( m_gpublocks == 0 ) throw std::runtime_error( "RamboSamplingKernelDevice: gpublocks must be > 0" );
if( m_gputhreads == 0 ) throw std::runtime_error( "RamboSamplingKernelDevice: gputhreads must be > 0" );
if( this->nevt() != m_rndmom.nevt() ) throw std::runtime_error( "RamboSamplingKernelDevice: nevt mismatch with rndmom" );
if( this->nevt() != m_momenta.nevt() ) throw std::runtime_error( "RamboSamplingKernelDevice: nevt mismatch with momenta" );
if( this->nevt() != m_weights.nevt() ) throw std::runtime_error( "RamboSamplingKernelDevice: nevt mismatch with weights" );
if( !m_rndmom.isOnDevice() ) throw std::runtime_error( "MasslessRamboSamplingKernelDevice: rndmom must be a device array" );
if( !m_momenta.isOnDevice() ) throw std::runtime_error( "MasslessRamboSamplingKernelDevice: momenta must be a device array" );
if( !m_weights.isOnDevice() ) throw std::runtime_error( "MasslessRamboSamplingKernelDevice: weights must be a device array" );
if( m_gpublocks == 0 ) throw std::runtime_error( "MasslessRamboSamplingKernelDevice: gpublocks must be > 0" );
if( m_gputhreads == 0 ) throw std::runtime_error( "MasslessRamboSamplingKernelDevice: gputhreads must be > 0" );
if( this->nevt() != m_rndmom.nevt() ) throw std::runtime_error( "MasslessRamboSamplingKernelDevice: nevt mismatch with rndmom" );
if( this->nevt() != m_momenta.nevt() ) throw std::runtime_error( "MasslessRamboSamplingKernelDevice: nevt mismatch with momenta" );
if( this->nevt() != m_weights.nevt() ) throw std::runtime_error( "MasslessRamboSamplingKernelDevice: nevt mismatch with weights" );
// Sanity checks for memory access (momenta buffer)
constexpr int neppM = MemoryAccessMomenta::neppM; // AOSOA layout
static_assert( ispoweroftwo( neppM ), "neppM is not a power of 2" );
if( m_gputhreads % neppM != 0 )
{
std::ostringstream sstr;
sstr << "RamboSamplingKernelHost: gputhreads should be a multiple of neppM=" << neppM;
sstr << "MasslessRamboSamplingKernelHost: gputhreads should be a multiple of neppM=" << neppM;
throw std::runtime_error( sstr.str() );
}
// Sanity checks for memory access (random number buffer)
Expand All @@ -128,7 +190,7 @@ namespace mg5amcCpu
if( m_gputhreads % neppR != 0 )
{
std::ostringstream sstr;
sstr << "RamboSamplingKernelDevice: gputhreads should be a multiple of neppR=" << neppR;
sstr << "MasslessRamboSamplingKernelDevice: gputhreads should be a multiple of neppR=" << neppR;
throw std::runtime_error( sstr.str() );
}
}
Expand All @@ -141,7 +203,7 @@ namespace mg5amcCpu
getMomentaInitialDevice( const fptype energy,
fptype* momenta )
{
constexpr auto getMomentaInitial = ramboGetMomentaInitial<DeviceAccessMomenta>;
constexpr auto getMomentaInitial = massless_rambo::ramboGetMomentaInitial<DeviceAccessMomenta>;
return getMomentaInitial( energy, momenta );
}
#endif
Expand All @@ -150,7 +212,7 @@ namespace mg5amcCpu

#ifdef MGONGPUCPP_GPUIMPL
void
RamboSamplingKernelDevice::getMomentaInitial()
MasslessRamboSamplingKernelDevice::getMomentaInitial()
{
gpuLaunchKernel( getMomentaInitialDevice, m_gpublocks, m_gputhreads, m_energy, m_momenta.data() );
}
Expand All @@ -165,7 +227,7 @@ namespace mg5amcCpu
fptype* momenta,
fptype* wgts )
{
constexpr auto getMomentaFinal = ramboGetMomentaFinal<DeviceAccessRandomNumbers, DeviceAccessMomenta, DeviceAccessWeights>;
constexpr auto getMomentaFinal = massless_rambo::ramboGetMomentaFinal<DeviceAccessRandomNumbers, DeviceAccessMomenta, DeviceAccessWeights>;
return getMomentaFinal( energy, rndmom, momenta, wgts );
}
#endif
Expand All @@ -174,7 +236,7 @@ namespace mg5amcCpu

#ifdef MGONGPUCPP_GPUIMPL
void
RamboSamplingKernelDevice::getMomentaFinal()
MasslessRamboSamplingKernelDevice::getMomentaFinal()
{
gpuLaunchKernel( getMomentaFinalDevice, m_gpublocks, m_gputhreads, m_energy, m_rndmom.data(), m_momenta.data(), m_weights.data() );
}
Expand Down
50 changes: 44 additions & 6 deletions madgraph/iolibs/template_files/madmatrix/RamboSamplingKernels.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,8 @@

#include "MemoryBuffers.h"

#include <vector>

#ifdef MGONGPUCPP_GPUIMPL
namespace mg5amcGpu
#else
Expand Down Expand Up @@ -68,19 +70,19 @@ namespace mg5amcCpu
//--------------------------------------------------------------------------

// A class encapsulating RAMBO phase space sampling on a CPU host
class RamboSamplingKernelHost final : public SamplingKernelBase, public NumberOfEvents
class MasslessRamboSamplingKernelHost final : public SamplingKernelBase, public NumberOfEvents
{
public:

// Constructor from existing input and output buffers
RamboSamplingKernelHost( const fptype energy, // input: energy
MasslessRamboSamplingKernelHost( const fptype energy, // input: energy
const BufferRndNumMomenta& rndmom, // input: random numbers in [0,1]
BufferMomenta& momenta, // output: momenta
BufferWeights& weights, // output: weights
const size_t nevt );

// Destructor
virtual ~RamboSamplingKernelHost() {}
virtual ~MasslessRamboSamplingKernelHost() {}

// Get momenta of initial state particles
void getMomentaInitial() override final;
Expand All @@ -94,22 +96,58 @@ namespace mg5amcCpu

//--------------------------------------------------------------------------

// Compability port mirroring the massless momenta implementation
// For now own RNG internally (or keep to match Fortran) TODO
// rndmom just for interface (maybe delete later) TODO
class RamboSamplingKernelHost final : public SamplingKernelBase, public NumberOfEvents
{
public:

RamboSamplingKernelHost( const fptype energy, // input: energy
const BufferRndNumMomenta& rndmom, // input: random [0,1] UNUSED
const std::vector<fptype>& masses, // input: external-leg masses
const int ninitial, // input: #n initial-state particles
const size_t nevt, // input: #n events
BufferMomenta& momenta, // output: momenta
BufferWeights& weights); // output: weights

virtual ~RamboSamplingKernelHost() {}

// No-op, kept to match the massless
void getMomentaInitial() override final;

// All the magic here
void getMomentaFinal() override final;

bool isOnDevice() const override final { return false; }

private:

// The EXTERNAL masses
std::vector<double> m_masses;

// The number of inital particles
const int m_ninitial;
};

//--------------------------------------------------------------------------

#ifdef MGONGPUCPP_GPUIMPL
// A class encapsulating RAMBO phase space sampling on a GPU device
class RamboSamplingKernelDevice final : public SamplingKernelBase, public NumberOfEvents
class MasslessRamboSamplingKernelDevice final : public SamplingKernelBase, public NumberOfEvents
{
public:

// Constructor from existing input and output buffers
RamboSamplingKernelDevice( const fptype energy, // input: energy
MasslessRamboSamplingKernelDevice( const fptype energy, // input: energy
const BufferRndNumMomenta& rndmom, // input: random numbers in [0,1]
BufferMomenta& momenta, // output: momenta
BufferWeights& weights, // output: weights
const size_t gpublocks,
const size_t gputhreads );

// Destructor
virtual ~RamboSamplingKernelDevice() {}
virtual ~MasslessRamboSamplingKernelDevice() {}

// Get momenta of initial state particles
void getMomentaInitial() override final;
Expand Down
Loading
Loading