Commit df748f60 authored by Simone Pagan Griso's avatar Simone Pagan Griso
Browse files

separate original file from repository

parent c76789b4
......@@ -57,6 +57,7 @@ atlas_add_library( Pythia8_iLib
src/UserHooks/WZPolarization.cxx
src/UserHooks/mHatReweight.cxx
src/UserHooks/DecayToSUEP.cxx
src/UserHooks/suep_shower.cxx
src/UserHooks/main31.cxx
src/UserHooks/EnhanceSplittings.cxx
src/UserResonances/ResonanceExcitedCI.cxx
......
......@@ -12,6 +12,8 @@
#include <stdexcept>
#include <iostream>
#include "suep_shower.h"
namespace Pythia8{
class DecayToSUEP;
}
......@@ -86,220 +88,6 @@ namespace Pythia8{
};
//******************************
//** Auxiliary class for SUEP generation
//******************************
class tolerance {
public:
tolerance(double eps) :
_eps(eps) {
}
bool operator()(double a, double b) {
return (fabs(b - a) <= _eps);
}
private:
double _eps;
};
class Suep_shower
{
public:
// Constructor
Suep_shower(double mass, double temperature, double energy);
/** Generate a shower event, in the rest frame of the showe. */
std::vector< std::vector <double> > generate_shower();
protected:
/** Mass of the dark mesons to be generated in the shower */
double m;
/** Temperature parameter */
double Temp;
/** Total energy of the system */
double Etot;
protected:
/** Maxwell-Boltzman distribution, slightly massaged */
double f(double p);
/** Derivative of Maxwell-Boltzmann */
double fp(double p);
/** Test function to be solved for p_plus,p_minus */
double test_fun(double p);
/** generate one random 4 vector from the thermal distribution */
std::vector<double> generateFourVector();
/** auxiliary function which computes the total energy difference as a function of the
* momentum vectors and a scale factor "a" to balance energy */
double reballance_func(double a, const std::vector< std::vector <double> > &event);
private:
/** mass/Temperature ratio */
double A;
/** convenience function of mass/Temperature ratio */
double p_m;
/** solutions for shower generation. See paper for details. */
double p_plus, p_minus;
/** solutions in Lambda for shower generation */
double lambda_plus,lambda_minus;
/** solutions for shower generation. See paper for details. */
double q_plus, q_minus, q_m;
};
using namespace boost::math::tools; // For bracket_and_solve_root.
// constructor
Suep_shower::Suep_shower(double mass, double temperature, double energy) {
m = mass;
Temp=temperature;
Etot=energy;
A=m/Temp;
p_m=std::sqrt(2/(A*A)*(1+std::sqrt(1+A*A)));
double pmax=std::sqrt(2+2*std::sqrt(1+A*A))/A; // compute the location of the maximum, to split the range
tolerance tol = 0.00001;
p_plus = (bisect(boost::bind(&Suep_shower::test_fun, this, _1),pmax,50.0, tol)).first; // first root
p_minus = (bisect(boost::bind(&Suep_shower::test_fun, this, _1), 0.0,pmax, tol)).first; // second root
lambda_plus = - f(p_plus)/fp(p_plus);
lambda_minus = f(p_minus)/fp(p_minus);
q_plus = lambda_plus / (p_plus - p_minus);
q_minus = lambda_minus / (p_plus - p_minus);
q_m = 1- (q_plus + q_minus);
}
// maxwell-boltzman distribution, slightly massaged
double Suep_shower::f(double p){
return p*p*exp(-A*p*p/(1+std::sqrt(1+p*p)));
}
// derivative of maxwell-boltzmann
double Suep_shower::fp(double p){
return exp(-A*p*p/(1+std::sqrt(1+p*p)))*p*(2-A*p*p/std::sqrt(1+p*p));
}
// test function to be solved for p_plus, p_minus
double Suep_shower::test_fun(double p){
return log(f(p)/f(p_m))+1.0;
}
// generate one random 4 vector from the thermal distribution
vector<double> Suep_shower::generateFourVector(){
vector<double> 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){
U = ((double) rand() / RAND_MAX);
V = ((double) rand() / RAND_MAX);
if(U < q_m){
Y=U/q_m;
X=( 1 - Y )*( p_minus + lambda_minus )+Y*( p_plus - lambda_plus );
if(V < f(X) / f(p_m) && X>0){
break;
}
}
else{if(U < q_m + q_plus){
E = -log((U-q_m)/q_plus);
X = p_plus - lambda_plus*(1-E);
if(V<exp(E)*f(X)/f(p_m) && X>0){
break;
}
}
else{
E = - log((U-(q_m+q_plus))/q_minus);
X = p_minus + lambda_minus * (1 - E);
if(V < exp(E)*f(X)/f(p_m) && X>0){
break;
}
}
}
}
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);
// compose the 4 vector
en = std::sqrt(p*p+(this->m)*(this->m));
fourvec.push_back(en);
fourvec.push_back(p*cos(phi)*sin(theta));
fourvec.push_back(p*sin(phi)*sin(theta));
fourvec.push_back(p*cos(theta));
return fourvec;
}
// auxiliary function which computes the total energy difference as a function of the momentum vectors and a scale factor "a"
// to ballance energy, we solve for "a" by demanding that this function vanishes
// By rescaling the momentum rather than the energy, I avoid having to do these annoying rotations from the previous version
double Suep_shower::reballance_func(double a, const vector< vector <double> > &event){
double result =0.0;
double p2;
for(unsigned n = 0; n<event.size();n++){
p2 = event[n][1]*event[n][1] + event[n][2]*event[n][2] + event[n][3]*event[n][3];
result += std::sqrt(a*a*p2 + (this->m)* (this->m));
}
return result - (this->Etot);
}
// generate a shower event, in the rest frame of the shower
vector< vector <double> > Suep_shower::generate_shower(){
vector<vector<double> > event;
double sum_E = 0.0;
// fill up event record
while(sum_E<(this->Etot)){
event.push_back(this->generateFourVector());
sum_E += (event.back()).at(0);
}
// 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;
}
}
// finally, ballance the total energy, without destroying momentum conservation
tolerance tol = 0.00001;
double p_scale;
p_scale = (bisect(boost::bind(&Suep_shower::reballance_func, this, _1, event),0.0,2.0, tol)).first;
for(int n=0;n<len;n++){
event[n][1] = p_scale*event[n][1];
event[n][2] = p_scale*event[n][2];
event[n][3] = p_scale*event[n][3];
// force the energy with the on-shell condition
event[n][0] = std::sqrt(event[n][1]*event[n][1] + event[n][2]*event[n][2] + event[n][3]*event[n][3] + (this->m)*(this->m));
}
return event;
}
//******************************
//** Implementation of DecayToSUEP UserHook
......@@ -377,12 +165,14 @@ namespace Pythia8{
} // loop over particles in the event
#ifdef SUEP_DEBUG
if (not particleFound) {
std::cout << "[SUEP_DEBUG] " << "Particle " << m_pdgId(settingsPtr) << " not found. Nothing to decay to SUEP for this event." << std::endl;
std::cout << "[DecayToSUEP] " << "Particle " << m_pdgId(settingsPtr) << " not found. Nothing to decay to SUEP for this event." << std::endl;
} else {
#ifdef SUEP_DEBUG
std::cout << "[SEUP_DEBUG] " << "All Done for this event." << std::endl;
#endif
}
#ifdef SUEP_DEBUG
std::cout << "[SUEP_DEBUG] Printing event after adding SUEP:" << std::endl;
for (int ii=0; ii < process.size(); ++ii) {
std::cout << "[SUEP_DEBUG] " << ii << ": id=" << process[ii].id() << ", Final=" << process[ii].isFinal() << ", mayDecay=" << process[ii].mayDecay() << ", Status=" << process[ii].status() << ", daughter1=" << process[ii].daughter1() << ", daughter2=" << process[ii].daughter2() << std::endl;
......
/*
* This file is taken from the available public code at:
* https://gitlab.com/simonknapen/suep_generator
* by Simon Knapen.
*/
#include "suep_shower.h"
using namespace boost::math::tools; // For bracket_and_solve_root.
using namespace std;
// constructor
Suep_shower::Suep_shower(double mass, double temperature, double energy) {
m = mass;
Temp=temperature;
Etot=energy;
A=m/Temp;
p_m=std::sqrt(2/(A*A)*(1+std::sqrt(1+A*A)));
double pmax=std::sqrt(2+2*std::sqrt(1+A*A))/A; // compute the location of the maximum, to split the range
tolerance tol = 0.00001;
p_plus = (bisect(boost::bind(&Suep_shower::test_fun, this, _1),pmax,50.0, tol)).first; // first root
p_minus = (bisect(boost::bind(&Suep_shower::test_fun, this, _1), 0.0,pmax, tol)).first; // second root
lambda_plus = - f(p_plus)/fp(p_plus);
lambda_minus = f(p_minus)/fp(p_minus);
q_plus = lambda_plus / (p_plus - p_minus);
q_minus = lambda_minus / (p_plus - p_minus);
q_m = 1- (q_plus + q_minus);
}
// maxwell-boltzman distribution, slightly massaged
double Suep_shower::f(double p){
return p*p*exp(-A*p*p/(1+std::sqrt(1+p*p)));
}
// derivative of maxwell-boltzmann
double Suep_shower::fp(double p){
return exp(-A*p*p/(1+std::sqrt(1+p*p)))*p*(2-A*p*p/std::sqrt(1+p*p));
}
// test function to be solved for p_plus, p_minus
double Suep_shower::test_fun(double p){
return log(f(p)/f(p_m))+1.0;
}
// generate one random 4 vector from the thermal distribution
vector<double> Suep_shower::generateFourVector(){
vector<double> 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){
U = ((double) rand() / RAND_MAX);
V = ((double) rand() / RAND_MAX);
if(U < q_m){
Y=U/q_m;
X=( 1 - Y )*( p_minus + lambda_minus )+Y*( p_plus - lambda_plus );
if(V < f(X) / f(p_m) && X>0){
break;
}
}
else{if(U < q_m + q_plus){
E = -log((U-q_m)/q_plus);
X = p_plus - lambda_plus*(1-E);
if(V<exp(E)*f(X)/f(p_m) && X>0){
break;
}
}
else{
E = - log((U-(q_m+q_plus))/q_minus);
X = p_minus + lambda_minus * (1 - E);
if(V < exp(E)*f(X)/f(p_m) && X>0){
break;
}
}
}
}
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);
// compose the 4 vector
en = std::sqrt(p*p+(this->m)*(this->m));
fourvec.push_back(en);
fourvec.push_back(p*cos(phi)*sin(theta));
fourvec.push_back(p*sin(phi)*sin(theta));
fourvec.push_back(p*cos(theta));
return fourvec;
}
// auxiliary function which computes the total energy difference as a function of the momentum vectors and a scale factor "a"
// to ballance energy, we solve for "a" by demanding that this function vanishes
// By rescaling the momentum rather than the energy, I avoid having to do these annoying rotations from the previous version
double Suep_shower::reballance_func(double a, const vector< vector <double> > &event){
double result =0.0;
double p2;
for(unsigned n = 0; n<event.size();n++){
p2 = event[n][1]*event[n][1] + event[n][2]*event[n][2] + event[n][3]*event[n][3];
result += std::sqrt(a*a*p2 + (this->m)* (this->m));
}
return result - (this->Etot);
}
// generate a shower event, in the rest frame of the shower
vector< vector <double> > Suep_shower::generate_shower(){
vector<vector<double> > event;
double sum_E = 0.0;
// fill up event record
while(sum_E<(this->Etot)){
event.push_back(this->generateFourVector());
sum_E += (event.back()).at(0);
}
// 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;
}
}
// finally, ballance the total energy, without destroying momentum conservation
tolerance tol = 0.00001;
double p_scale;
p_scale = (bisect(boost::bind(&Suep_shower::reballance_func, this, _1, event),0.0,2.0, tol)).first;
for(int n=0;n<len;n++){
event[n][1] = p_scale*event[n][1];
event[n][2] = p_scale*event[n][2];
event[n][3] = p_scale*event[n][3];
// force the energy with the on-shell condition
event[n][0] = std::sqrt(event[n][1]*event[n][1] + event[n][2]*event[n][2] + event[n][3]*event[n][3] + (this->m)*(this->m));
}
return event;
}
/*
* This file is taken from the available public code at:
* https://gitlab.com/simonknapen/suep_generator
* by Simon Knapen.
*/
#ifndef SUEPSHOWER_h
#define SUEPSHOWER_h
#include <vector>
#include <math.h>
#include <boost/math/tools/roots.hpp>
#include <boost/bind.hpp>
//** Auxiliary class for tolerance checks */
class tolerance {
public:
tolerance(double eps) :
_eps(eps) {
}
bool operator()(double a, double b) {
return (fabs(b - a) <= _eps);
}
private:
double _eps;
};
/** Auxiliary class for SUEP generation.
* Details on models available on arXiv:1612.00850.
*/
class Suep_shower
{
public:
// Constructor
Suep_shower(double mass, double temperature, double energy);
/** Generate a shower event, in the rest frame of the showe. */
std::vector< std::vector <double> > generate_shower();
protected:
/** Mass of the dark mesons to be generated in the shower */
double m;
/** Temperature parameter */
double Temp;
/** Total energy of the system */
double Etot;
protected:
/** Maxwell-Boltzman distribution, slightly massaged */
double f(double p);
/** Derivative of Maxwell-Boltzmann */
double fp(double p);
/** Test function to be solved for p_plus,p_minus */
double test_fun(double p);
/** generate one random 4 vector from the thermal distribution */
std::vector<double> generateFourVector();
/** auxiliary function which computes the total energy difference as a function of the
* momentum vectors and a scale factor "a" to balance energy */
double reballance_func(double a, const std::vector< std::vector <double> > &event);
private:
/** mass/Temperature ratio */
double A;
/** convenience function of mass/Temperature ratio */
double p_m;
/** solutions for shower generation. See paper for details. */
double p_plus, p_minus;
/** solutions in Lambda for shower generation */
double lambda_plus,lambda_minus;
/** solutions for shower generation. See paper for details. */
double q_plus, q_minus, q_m;
};
#endif
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