Commit 5ffde251 authored by Simone Pagan Griso's avatar Simone Pagan Griso
Browse files

move to Pythia8 random number generator

parent 91dbb383
......@@ -115,7 +115,7 @@ namespace Pythia8{
particleFound=true;
//setup SUEP shower
static Suep_shower suep_shower(m_darkMesonMass(settingsPtr), m_darkTemperature(settingsPtr), m_mass(settingsPtr));
static Suep_shower suep_shower(m_darkMesonMass(settingsPtr), m_darkTemperature(settingsPtr), m_mass(settingsPtr), rndmPtr);
#ifdef SUEP_DEBUG
std::cout << "[SUEP_DEBUG] " << "Particle (pdgId=" << m_pdgId(settingsPtr) << ", isFinal=True) found. Decaying to SUEP now." << std::endl;
......
......@@ -10,10 +10,11 @@ using namespace std;
using namespace Pythia8;
// constructor
Suep_shower::Suep_shower(double mass, double temperature, double energy) {
Suep_shower::Suep_shower(double mass, double temperature, double energy, Pythia8::Rndm *rndm) {
m = mass;
Temp=temperature;
Etot=energy;
rndmEngine = rndm;
A=m/Temp;
p_m=std::sqrt(2/(A*A)*(1+std::sqrt(1+A*A)));
......@@ -56,8 +57,13 @@ using namespace Pythia8;
double U(0.0), V(0.0), X(0.0), Y(0.0), E(0.0);
int i=0;
while(i<100){
U = ((double) rand() / RAND_MAX);
V = ((double) rand() / RAND_MAX);
if (rndmEngine) {
U = rndmEngine->flat();
V = rndmEngine->flat();
} else {
U = ((double) rand() / RAND_MAX);
V = ((double) rand() / RAND_MAX);
}
if(U < q_m){
Y=U/q_m;
......@@ -85,8 +91,13 @@ using namespace Pythia8;
p=X*(this->m); // X is the dimensionless momentum, p/m
// now do the angles
phi = 2.0*M_PI*((double) rand() / RAND_MAX);
theta = acos(2.0*((double) rand() / RAND_MAX)-1.0);
if (rndmEngine) {
phi = 2.0*M_PI*(rndmEngine->flat());
theta = acos(2.0*rndmEngine->flat()-1.0);
} else {
phi = 2.0*M_PI*((double) rand() / RAND_MAX);
theta = acos(2.0*((double) rand() / RAND_MAX)-1.0);
}
// compose the 4 vector
en = std::sqrt(p*p+(this->m)*(this->m));
......
......@@ -11,6 +11,7 @@
#include <boost/math/tools/roots.hpp>
#include <boost/bind.hpp>
#include "Pythia8/PhaseSpace.h"
#include "Pythia8/Basics.h"
//** Auxiliary class for tolerance checks */
class tolerance {
......@@ -32,13 +33,21 @@ class Suep_shower
{
public:
// Constructor
Suep_shower(double mass, double temperature, double energy);
/** Constructor.
* @param mass Mass of the dark meson
* @param temperature model parameter
* @param energy total energy of decaying system
* @param rndm random number generator, if any (if not provided will use the rand() function).
*/
Suep_shower(double mass, double temperature, double energy, Pythia8::Rndm *rndm = 0);
/** Generate a shower event, in the rest frame of the showe. */
std::vector< Pythia8::Vec4 > generate_shower();
protected:
/** Random number generator, if not provided will use rand() */
Pythia8::Rndm *rndmEngine;
/** Mass of the dark mesons to be generated in the shower */
double m;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment