Suep_shower Class Reference

Auxiliary class for SUEP generation. More...

#include <suep_shower.h>

Collaboration diagram for Suep_shower:

Public Member Functions
	Suep_shower (double mass, double temperature, double energy, Pythia8::Rndm *rndm=0)
	Constructor.
std::vector< Pythia8::Vec4 >	generate_shower ()
	Generate a shower event, in the rest frame of the showe.

Protected Member Functions
double	f (double p)
	Maxwell-Boltzman distribution, slightly massaged.
double	fp (double p)
	Derivative of Maxwell-Boltzmann.
double	test_fun (double p)
	Test function to be solved for p_plus,p_minus.
Pythia8::Vec4	generateFourVector ()
	generate one random 4 vector from the thermal distribution
double	reballance_func (double a, const std::vector< Pythia8::Vec4 > &event)
	auxiliary function which computes the total energy difference as a function of the momentum vectors and a scale factor "a" to balance energy

Protected Attributes
Pythia8::Rndm *	m_rndmEngine
	Random number generator, if not provided will use rand().
double	m_m
	Mass of the dark mesons to be generated in the shower.
double	m_Temp
	Temperature parameter.
double	m_Etot
	Total energy of the system.

Private Attributes
double	m_A
	mass/Temperature ratio
double	m_p_m
	convenience function of mass/Temperature ratio
double	m_p_plus
	solutions for shower generation.
double	m_p_minus
double	m_lambda_plus
	solutions in Lambda for shower generation
double	m_lambda_minus
double	m_q_plus
	solutions for shower generation.
double	m_q_minus
double	m_q_m

Detailed Description

Auxiliary class for SUEP generation.

Details on models available on arXiv:1612.00850.

Definition at line 30 of file suep_shower.h.

Constructor & Destructor Documentation

Suep_shower()

Suep_shower::Suep_shower	(	double	mass,
		double	temperature,
		double	energy,
		Pythia8::Rndm *	rndm = 0 )

Constructor.

Parameters

mass	Mass of the dark meson
temperature	model parameter
energy	total energy of decaying system
rndm	random number generator, if any (if not provided will use the rand() function).

Definition at line 15 of file suep_shower.cxx.

                                                                                          {
  m_m = mass;
  m_Temp=temperature;
  m_Etot=energy;
  m_rndmEngine = rndm;
  
  m_A=m_m/m_Temp;
  m_p_m=std::sqrt(2/(m_A*m_A)*(1+std::sqrt(1+m_A*m_A)));
  
  double pmax=std::sqrt(2+2*std::sqrt(1+m_A*m_A))/m_A; // compute the location of the maximum, to split the range 
  
  tolerance tol = 0.00001;
  m_p_plus = (bisect(std::bind(&Suep_shower::test_fun, this, std::placeholders::_1),pmax,50.0, tol)).first; // first root
  m_p_minus = (bisect(std::bind(&Suep_shower::test_fun, this, std::placeholders::_1), 0.0,pmax, tol)).first; // second root
  m_lambda_plus = - f(m_p_plus)/fp(m_p_plus);
  m_lambda_minus = f(m_p_minus)/fp(m_p_minus);
  m_q_plus = m_lambda_plus / (m_p_plus - m_p_minus);
  m_q_minus = m_lambda_minus / (m_p_plus - m_p_minus);
  m_q_m = 1- (m_q_plus + m_q_minus);
  
}

Member Function Documentation

f()

double Suep_shower::f ( double p )

protected

Maxwell-Boltzman distribution, slightly massaged.

Definition at line 38 of file suep_shower.cxx.

                              {
  return p*p*exp(-m_A*p*p/(1+std::sqrt(1+p*p)));
}

fp()

double Suep_shower::fp ( double p )

protected

Derivative of Maxwell-Boltzmann.

Definition at line 43 of file suep_shower.cxx.

                               {
  return exp(-m_A*p*p/(1+std::sqrt(1+p*p)))*p*(2-m_A*p*p/std::sqrt(1+p*p));
}

generate_shower()

vector< Vec4 > Suep_shower::generate_shower

(

)

Generate a shower event, in the rest frame of the showe.

Definition at line 128 of file suep_shower.cxx.

                                           {
  
  vector< Vec4 > event;
  double sum_E = 0.0;
  
  // fill up event record
  while(sum_E<(this->m_Etot)){
    event.push_back(this->generateFourVector());
    sum_E += (event.back()).e();
  }
  
  // reballance momenta
  int len = event.size();
  double sum_p, correction;
  for(int i = 1;i<4;i++){ // loop over 3 carthesian directions
    
    sum_p = 0.0;
    for(int n=0;n<len;n++){
      sum_p+=event[n][i];
    }
    correction=-1.0*sum_p/len;
    
    for(int n=0;n<len;n++){
      event[n][i] += correction;
    } 
  }
 
  //Shield against an exception in the calculation of "p_scale" further down. If it fails, abort the event.
  if(Suep_shower::reballance_func(2.0,event)<0.0){
    // failed to balance energy. 
    event.clear();
    return event;
  } 
  
  // finally, ballance the total energy, without destroying momentum conservation
  tolerance tol = 0.00001;
  double p_scale;
  try {
    p_scale = (bisect(std::bind(&Suep_shower::reballance_func, this, std::placeholders::_1, event),0.0,2.0, tol)).first;
  } catch (std::exception &e) {
    //in some rare circumstances, this balancing might fail
    std::cout << "[SUEP_SHOWER] WARNING: Failed to rebalance the following event; Printing for debug and throw exception" << std::endl;
    std::cout << e.what() << std::endl;
    std::cout << "N. Particle, px, py, pz, E" << std::endl;
    for (size_t jj=0; jj < event.size(); jj++) {
      //print out event
      std::cout << jj << ": " << event[jj].px() << ", " << event[jj].py() << ", " << event[jj].pz() << ", " << event[jj].e() << std::endl;
    }
    throw e;
  }
  
  for(int n=0;n<len;n++){
    event[n].px(p_scale*event[n].px());
    event[n].py(p_scale*event[n].py());
    event[n].pz(p_scale*event[n].pz());
    // force the energy with the on-shell condition
    event[n].e(std::sqrt(event[n].px()*event[n].px() + event[n].py()*event[n].py() + event[n].pz()*event[n].pz() + (this->m_m)*(this->m_m)));
  }
  
  return event;    
}

generateFourVector()

Vec4 Suep_shower::generateFourVector ( )

protected

generate one random 4 vector from the thermal distribution

Definition at line 53 of file suep_shower.cxx.

                                    {
  
  Vec4 fourvec;
  double en, phi, theta, p;//kinematic variables of the 4 vector
  
  // first do momentum, following arxiv:1305.5226
  double U(0.0), V(0.0), X(0.0), Y(0.0), E(0.0);
  int i=0;      
  while(i<100){
    if (m_rndmEngine) {
      U = m_rndmEngine->flat();
      V = m_rndmEngine->flat();
    } else {
      U = ((double) rand() / RAND_MAX);
      V = ((double) rand() / RAND_MAX);
    }
    
    if(U < m_q_m){
      Y=U/m_q_m;
      X=( 1 - Y )*( m_p_minus + m_lambda_minus )+Y*( m_p_plus - m_lambda_plus );
      if(V < f(X) / f(m_p_m) && X>0){
    break;
      }
    }
    else{if(U < m_q_m + m_q_plus){
    E = -log((U-m_q_m)/m_q_plus);
    X = m_p_plus - m_lambda_plus*(1-E);
    if(V<exp(E)*f(X)/f(m_p_m) && X>0){
      break;
    }
      }
      else{
    E = - log((U-(m_q_m+m_q_plus))/m_q_minus);
    X = m_p_minus + m_lambda_minus * (1 - E);
    if(V < exp(E)*f(X)/f(m_p_m) && X>0){
      break;
    }
      }
    }
  }
  p=X*(this->m_m); // X is the dimensionless momentum, p/m
  
  // now do the angles
  if (m_rndmEngine) {
    phi = 2.0*M_PI*(m_rndmEngine->flat());
    theta = acos(2.0*m_rndmEngine->flat()-1.0);
  } else {
    //coverity[DC.WEAK_CRYPTO]
    phi = 2.0*M_PI*((double) rand() / RAND_MAX);
    //coverity[DC.WEAK_CRYPTO]
    theta = acos(2.0*((double) rand() / RAND_MAX)-1.0);
  }
  
  // compose the 4 vector
  en = std::sqrt(p*p+(this->m_m)*(this->m_m));
  //coverity[COPY_PASTE_ERROR]
  fourvec.p(p*cos(phi)*sin(theta), p*sin(phi)*sin(theta), p*cos(theta), en);
  
  return fourvec; 
}

reballance_func()

double Suep_shower::reballance_func ( double a, const std::vector< Pythia8::Vec4 > & event )

protected

auxiliary function which computes the total energy difference as a function of the momentum vectors and a scale factor "a" to balance energy

Definition at line 117 of file suep_shower.cxx.

                                                                        {
  double result =0.0;
  double p2;
  for(unsigned n = 0; n<event.size();n++){
    p2 = event[n].px()*event[n].px() + event[n].py()*event[n].py() + event[n].pz()*event[n].pz();
    result += std::sqrt(a*a*p2 + (this->m_m)* (this->m_m));
  }
  return result - (this->m_Etot);
}

test_fun()

double Suep_shower::test_fun ( double p )

protected

Test function to be solved for p_plus,p_minus.

Definition at line 48 of file suep_shower.cxx.

                                     {
  return log(f(p)/f(m_p_m))+1.0;
}

Member Data Documentation

m_A

double Suep_shower::m_A

private

mass/Temperature ratio

Definition at line 79 of file suep_shower.h.

m_Etot

double Suep_shower::m_Etot

protected

Total energy of the system.

Definition at line 57 of file suep_shower.h.

m_lambda_minus

double Suep_shower::m_lambda_minus

private

Definition at line 85 of file suep_shower.h.

m_lambda_plus

double Suep_shower::m_lambda_plus

private

solutions in Lambda for shower generation

Definition at line 85 of file suep_shower.h.

m_m

double Suep_shower::m_m

protected

Mass of the dark mesons to be generated in the shower.

Definition at line 51 of file suep_shower.h.

m_p_m

double Suep_shower::m_p_m

private

convenience function of mass/Temperature ratio

Definition at line 81 of file suep_shower.h.

m_p_minus

double Suep_shower::m_p_minus

private

Definition at line 83 of file suep_shower.h.

m_p_plus

double Suep_shower::m_p_plus

private

solutions for shower generation.

See paper for details.

Definition at line 83 of file suep_shower.h.

m_q_m

double Suep_shower::m_q_m

private

Definition at line 87 of file suep_shower.h.

m_q_minus

double Suep_shower::m_q_minus

private

Definition at line 87 of file suep_shower.h.

m_q_plus

double Suep_shower::m_q_plus

private

solutions for shower generation.

See paper for details.

Definition at line 87 of file suep_shower.h.

m_rndmEngine

Pythia8::Rndm* Suep_shower::m_rndmEngine

protected

Random number generator, if not provided will use rand().

Definition at line 48 of file suep_shower.h.

m_Temp

double Suep_shower::m_Temp

protected

Temperature parameter.

Definition at line 54 of file suep_shower.h.

---

The documentation for this class was generated from the following files:

• suep_shower.h
• suep_shower.cxx