Page MenuHomeHEPForge
Files F7879443

No OneTemporary
Actions

View Options

diff --git a/CepGen/Core/Exception.h b/CepGen/Core/Exception.h
index 05b6447..f1aa770 100644
--- a/CepGen/Core/Exception.h
+++ b/CepGen/Core/Exception.h
@@ -1,134 +1,134 @@
#ifndef CepGen_Core_Exception_h
#define CepGen_Core_Exception_h
-#include <iomanip>
#include <sstream>
+#include <stdexcept>
#include "Logger.h"
#define PrintMessage( m ) \
if ( CepGen::Logger::get().level > CepGen::Logger::Nothing ) { CepGen::Exception( __PRETTY_FUNCTION__, m, CepGen::Verbatim ).dump( CepGen::Logger::get().outputStream ); }
#define Information( m ) \
if ( CepGen::Logger::get().level > CepGen::Logger::Nothing ) { CepGen::Exception( __PRETTY_FUNCTION__, m, CepGen::Information ).dump( CepGen::Logger::get().outputStream ); }
#define Debugging( m ) \
if ( CepGen::Logger::get().level >= CepGen::Logger::Debug ) { CepGen::Exception( __PRETTY_FUNCTION__, m, CepGen::DebugMessage ).dump( CepGen::Logger::get().outputStream ); }
#define DebuggingInsideLoop( m ) \
if ( CepGen::Logger::get().level >= CepGen::Logger::DebugInsideLoop ) { CepGen::Exception( __PRETTY_FUNCTION__, m, CepGen::DebugMessage ).dump( CepGen::Logger::get().outputStream ); }
#define InWarning( m ) \
if ( CepGen::Logger::get().level >= CepGen::Logger::Warning ) { CepGen::Exception( __PRETTY_FUNCTION__, m, CepGen::JustWarning ).dump( CepGen::Logger::get().outputStream ); }
#define InError( m ) \
if ( CepGen::Logger::get().level >= CepGen::Logger::Error ) { CepGen::Exception( __PRETTY_FUNCTION__, m, CepGen::ErrorMessage ).dump( CepGen::Logger::get().outputStream ); }
#define FatalError( m ) \
CepGen::Exception( __PRETTY_FUNCTION__, m, CepGen::FatalError ).dump( CepGen::Logger::get().outputStream );
namespace CepGen
{
/**
* \brief Enumeration of exception severities
* \author Laurent Forthomme <laurent.forthomme@cern.ch>
* \date 27 Mar 2015
*/
enum ExceptionType { Undefined=-1, Verbatim, Information, DebugMessage, JustWarning, ErrorMessage, FatalError };
/**
* \brief A simple exception handler
* \author Laurent Forthomme <laurent.forthomme@cern.ch>
* \date 24 Mar 2015
*/
class Exception : public std::runtime_error
{
public:
/// Initialize a new exception object
/// \param[in] from method invoking the exception
/// \param[in] desc brief description of the exception
/// \param[in] type exception type
/// \param[in] id exception code (useful for logging)
inline Exception( const char* from, std::string desc, ExceptionType type=Undefined, const int id=0 ) :
std::runtime_error( desc ), from_( from ), type_( type ), error_num_( id ) {}
/// Initialize a new exception object
/// \param[in] from method invoking the exception
/// \param[in] desc brief description of the exception
/// \param[in] type exception type
/// \param[in] id exception code (useful for logging)
inline Exception( const char* from, const char* desc, ExceptionType type=Undefined, const int id=0 ) :
std::runtime_error( desc ), from_( from ), type_( type ), error_num_( id ) {}
inline ~Exception() throw() {
if ( type() == FatalError ) exit(0);
// we stop this process' execution on fatal exception
}
/// Extract the origin of the exception
inline std::string from() const { return from_; }
/// Extract the exception code
inline int errorNumber() const { return error_num_; }
/// Extract the exception type
inline ExceptionType type() const { return type_; }
/// Extract a human-readable (and colourified) version of the exception type
inline std::string typeString() const {
switch ( type() ) {
case JustWarning: return "\033[34;1mJustWarning\033[0m";
case Information: return "\033[32;1mInfo\033[0m";
case DebugMessage: return "\033[33;1mDebug\033[0m";
case ErrorMessage: return "\033[31;1mError\033[0m";
case FatalError: return "\033[31;1mFatal\033[0m";
case Undefined: default: return "\33[7;1mUndefined\033[0m";
}
}
/// Dump the full exception information in a given output stream
/// \param[inout] os the output stream where the information is dumped
- inline void dump(std::ostream& os=std::cerr) const {
+ inline void dump( std::ostream& os = Logger::get().outputStream ) const {
if ( type() == Verbatim ) {
os << what() << std::endl;
return;
}
if ( type() == Information ) {
os << "[\033[32;1mInfo.\033[0m]\t" << what() << std::endl;
return;
}
if ( type() == DebugMessage ) {
os << "==================================== \033[33;1mDebug\033[0m ====================================" << std::endl
<< " From: " << from() << std::endl;
}
else {
os << "============================= Exception detected! =============================" << std::endl
<< " Class: " << typeString() << std::endl
<< " Raised by: " << from() << std::endl;
}
os << " Description: \t" << what() << std::endl;
if ( errorNumber() != 0 )
os << "-------------------------------------------------------------------------------" << std::endl
<< " Error #" << errorNumber() << std::endl;
os << "===============================================================================" << std::endl;
}
/// Extract a one-line summary of the exception
inline std::string OneLine() const {
std::ostringstream os;
os << "[" << type() << "] === " << from() << " === "
<< what();
return os.str();
}
private:
/// Origin of the exception
std::string from_;
/// Exception type
ExceptionType type_;
/// Integer exception number
int error_num_;
};
class Printer
{
public:
inline static Exception LogInfo( const char* name ) { return Exception( "", name, Information ); }
void operator<<( const char* text ) {
LogInfo( text ).dump( Logger::get().outputStream );
}
};
}
#endif
diff --git a/CepGen/Core/Generator.cpp b/CepGen/Core/Generator.cpp
index a62dca8..72105bc 100644
--- a/CepGen/Core/Generator.cpp
+++ b/CepGen/Core/Generator.cpp
@@ -1,115 +1,138 @@
#include "CepGen/Generator.h"
+#include "CepGen/Parameters.h"
+
+#include "CepGen/Core/Vegas.h"
+
+#include "CepGen/Event/Event.h"
#include <fstream>
namespace CepGen
{
Generator::Generator() :
parameters( std::unique_ptr<Parameters>( new Parameters ) ),
cross_section_( -1. ), cross_section_error_( -1. ), has_cross_section_( false )
{
if ( Logger::get().level > Logger::Nothing ) {
Debugging( "Generator initialized" );
try { printHeader(); } catch ( Exception& e ) { e.dump(); }
}
srand( time( 0 ) ); // Random number initialization
}
Generator::Generator( Parameters* ip ) :
parameters( ip )
{}
Generator::~Generator()
{
if ( parameters->generation.enabled && parameters->process() && parameters->process()->numGeneratedEvents()>0 ) {
Information( Form( "Mean generation time / event: %.3f ms", parameters->process()->totalGenerationTime()*1.e3/parameters->process()->numGeneratedEvents() ) );
}
}
+ size_t
+ Generator::numDimensions() const
+ {
+ if ( !parameters->process() ) return 0;
+ return parameters->process()->numDimensions( parameters->kinematics.mode );
+ }
+
void
Generator::clearRun()
{
vegas_.reset();
parameters->vegas.first_run = true;
has_cross_section_ = false; // force the recreation of the Vegas instance
cross_section_ = cross_section_error_ = -1.;
}
void
Generator::setParameters( Parameters& ip )
{
parameters = std::unique_ptr<Parameters>( new Parameters( ip ) ); // copy constructor
}
void
Generator::printHeader()
{
std::string tmp;
std::ostringstream os; os << std::endl;
std::ifstream hf( "README" );
if ( !hf.good() ) throw Exception( __PRETTY_FUNCTION__, "Failed to open README file", JustWarning );
while ( true ) {
if ( !hf.good() ) break;
getline( hf, tmp );
os << "\n " << tmp;
}
hf.close();
Information( os.str().c_str() );
}
+ double
+ Generator::computePoint( double* x )
+ {
+ prepareFunction();
+ double res = f( x, numDimensions(), (void*)parameters.get() );
+ std::ostringstream os;
+ for ( unsigned int i=0; i<numDimensions(); i++ ) { os << x[i] << " "; }
+ Debugging( Form( "Result for x[%zu] = ( %s):\n\t%10.6f", numDimensions(), os.str().c_str(), res ) );
+ return res;
+ }
+
void
Generator::computeXsection( double& xsec, double& err )
{
// first destroy and recreate the Vegas instance
if ( !vegas_ ) {
vegas_ = std::unique_ptr<Vegas>( new Vegas( numDimensions(), f, parameters.get() ) );
}
else if ( vegas_->dimensions() != numDimensions() ) {
vegas_.reset( new Vegas( numDimensions(), f, parameters.get() ) );
}
if ( Logger::get().level>=Logger::Debug ) {
std::ostringstream topo; topo << parameters->kinematics.mode;
Debugging( Form( "New Vegas instance created\n\t"
"Considered topology: %s case\n\t"
"Will proceed with %d-dimensional integration", topo.str().c_str(), numDimensions() ) );
}
Information( "Starting the computation of the process cross-section" );
try { prepareFunction(); } catch ( Exception& e ) { e.dump(); }
has_cross_section_ = ( vegas_->integrate( cross_section_, cross_section_error_ ) == 0 );
xsec = cross_section_;
err = cross_section_error_;
Information( Form( "Total cross section: %f +/- %f pb", xsec, err ) );
}
Event*
Generator::generateOneEvent()
{
bool good = false;
if ( !has_cross_section_ ) {
computeXsection( cross_section_, cross_section_error_ );
}
while ( !good ) { good = vegas_->generateOneEvent(); }
last_event = this->parameters->generation.last_event;
return last_event.get();
}
void
Generator::prepareFunction()
{
if ( !parameters->process() ) {
throw Exception( __PRETTY_FUNCTION__, "No process defined!", FatalError );
}
Kinematics kin = parameters->kinematics;
parameters->process()->addEventContent();
parameters->process()->setKinematics( kin );
Debugging( "Function prepared to be integrated!" );
}
}
diff --git a/CepGen/Core/Vegas.cpp b/CepGen/Core/Vegas.cpp
index ed9fdd8..5cf0782 100644
--- a/CepGen/Core/Vegas.cpp
+++ b/CepGen/Core/Vegas.cpp
@@ -1,314 +1,315 @@
#include "Vegas.h"
+#include "CepGen/Parameters.h"
#include <fstream>
namespace CepGen
{
Vegas::Vegas( const unsigned int dim, double f_( double*, size_t, void* ), Parameters* param ) :
vegas_bin_( 0 ), correc_( 0. ), correc2_( 0. ),
input_params_( param ),
grid_prepared_( false ), gen_prepared_( false ),
f_max2_( 0. ), f_max_diff_( 0. ), f_max_old_( 0. ), f_max_global_( 0. ),
function_( std::unique_ptr<gsl_monte_function>( new gsl_monte_function ) )
{
//--- function to be integrated
function_->f = f_;
function_->dim = dim;
function_->params = (void*)param;
num_converg_ = param->vegas.ncvg;
num_iter_ = param->vegas.itvg;
//--- initialise the random number generator
gsl_rng_env_setup();
rng_ = gsl_rng_alloc( gsl_rng_default );
unsigned long seed = ( param->vegas.seed > 0 )
? param->vegas.seed
: time( nullptr ); // seed with time
gsl_rng_set( rng_, seed );
Debugging( Form( "Number of integration dimensions: %d\n\t"
"Number of iterations: %d\n\t"
"Number of function calls: %d", dim, num_iter_, num_converg_ ) );
}
Vegas::~Vegas()
{
if ( rng_ ) gsl_rng_free( rng_ );
}
int
Vegas::integrate( double& result, double& abserr )
{
//--- prepare Vegas
gsl_monte_vegas_state* state = gsl_monte_vegas_alloc( function_->dim );
//--- integration bounds
std::vector<double> x_low( function_->dim, 0. ), x_up( function_->dim, 1. );
//--- launch Vegas
int veg_res;
//----- warmup (prepare the grid)
if ( !grid_prepared_ ) {
veg_res = gsl_monte_vegas_integrate( function_.get(), &x_low[0], &x_up[0], function_->dim, 10000, rng_, state, &result, &abserr );
grid_prepared_ = true;
}
//----- integration
for ( unsigned int i=0; i<num_iter_; i++ ) {
veg_res = gsl_monte_vegas_integrate( function_.get(), &x_low[0], &x_up[0], function_->dim, 0.2*num_converg_, rng_, state, &result, &abserr );
PrintMessage( Form( ">> Iteration %2d: average = %10.6f sigma = %10.6f chi2 = %4.3f", i+1, result, abserr, gsl_monte_vegas_chisq( state ) ) );
}
//--- clean Vegas
gsl_monte_vegas_free( state );
return veg_res;
}
void
Vegas::generate()
{
std::ofstream of;
std::string fn;
if ( !gen_prepared_ ) setGen();
Information( Form( "%d events will be generated", input_params_->generation.maxgen ) );
unsigned int i = 0;
while ( i < input_params_->generation.maxgen ) {
if ( generateOneEvent() ) i++;
}
Information( Form( "%d events generated", i ) );
}
bool
Vegas::generateOneEvent()
{
if ( !gen_prepared_ ) setGen();
const unsigned int ndim = function_->dim, max = pow( mbin_, ndim );
std::vector<double> x( ndim, 0. );
//--- correction cycles
if ( vegas_bin_ != 0 ) {
bool has_correction = false;
while ( !correctionCycle( x, has_correction ) ) {}
if ( has_correction ) return storeEvent( x );
}
double weight;
double y = -1.;
//--- normal generation cycle
//----- select a Vegas bin and reject if fmax is too small
do {
do {
// ...
vegas_bin_ = uniform() * max;
y = uniform() * f_max_global_;
nm_[vegas_bin_] += 1;
} while ( y > f_max_[vegas_bin_] );
// Select x values in this Vegas bin
int jj = vegas_bin_;
for ( unsigned int i=0; i<ndim; i++ ) {
int jjj = jj / mbin_;
n_[i] = jj - jjj * mbin_;
x[i] = ( uniform() + n_[i] ) / mbin_;
jj = jjj;
}
// Get weight for selected x value
weight = F( x );
} while ( y > weight );
if ( weight <= f_max_[vegas_bin_] ) vegas_bin_ = 0;
// Init correction cycle if weight is higher than fmax or ffmax
else if ( weight <= f_max_global_ ) {
f_max_old_ = f_max_[vegas_bin_];
f_max_[vegas_bin_] = weight;
f_max_diff_ = weight-f_max_old_;
correc_ = ( nm_[vegas_bin_] - 1. ) * f_max_diff_ / f_max_global_ - 1.;
}
else {
f_max_old_ = f_max_[vegas_bin_];
f_max_[vegas_bin_] = weight;
f_max_diff_ = weight-f_max_old_;
f_max_global_ = weight;
correc_ = ( nm_[vegas_bin_] - 1. ) * f_max_diff_ / f_max_global_ * weight / f_max_global_ - 1.;
}
Debugging( Form( "Correction applied: %f, Vegas bin = %d", correc_, vegas_bin_ ) );
// Return with an accepted event
if ( weight > 0. ) return storeEvent( x );
return false;
}
bool
Vegas::correctionCycle( std::vector<double>& x, bool& has_correction )
{
double weight;
const unsigned int ndim = function_->dim;
Debugging( Form( "Correction cycles are started.\n\t"
"j = %f"
"correc = %f"
"corre2 = %f", vegas_bin_, correc2_ ) );
if ( correc_ >= 1. ) correc_ -= 1.;
if ( uniform() < correc_ ) {
correc_ = -1.;
std::vector<double> xtmp( ndim, 0. );
// Select x values in Vegas bin
for ( unsigned int k=0; k<ndim; k++ ) {
xtmp[k] = ( uniform() + n_[k] ) * inv_mbin_;
}
// Compute weight for x value
weight = F( xtmp );
// Parameter for correction of correction
if ( weight > f_max_[vegas_bin_] ) {
if ( weight > f_max2_ ) f_max2_ = weight;
correc2_ -= 1.;
correc_ += 1.;
}
// Accept event
if ( weight >= f_max_diff_*uniform() + f_max_old_ ) { // FIXME!!!!
//Error("Accepting event!!!");
//return storeEvent(x);
x = xtmp;
has_correction = true;
return true;
}
return false;
}
// Correction if too big weight is found while correction
// (All your bases are belong to us...)
if ( f_max2_ > f_max_[vegas_bin_] ) {
f_max_old_ = f_max_[vegas_bin_];
f_max_[vegas_bin_] = f_max2_;
f_max_diff_ = f_max2_-f_max_old_;
if ( f_max2_ < f_max_global_ ) {
correc_ = ( nm_[vegas_bin_] - 1. ) * f_max_diff_ / f_max_global_ - correc2_;
}
else {
f_max_global_ = f_max2_;
correc_ = ( nm_[vegas_bin_] - 1. ) * f_max_diff_ / f_max_global_ * f_max2_ / f_max_global_ - correc2_;
}
correc2_ = 0.;
f_max2_ = 0.;
return false;
}
return true;
}
bool
Vegas::storeEvent( const std::vector<double>& x )
{
input_params_->setStorage( true );
F( x );
input_params_->generation.ngen += 1;
input_params_->setStorage( false );
if ( input_params_->generation.ngen % input_params_->generation.gen_print_every == 0 ) {
Debugging( Form( "Generated events: %d", input_params_->generation.ngen ) );
input_params_->generation.last_event->dump();
}
return true;
}
void
Vegas::setGen()
{
Information( Form( "Preparing the grid for the generation of unweighted events: %d points", input_params_->vegas.npoints ) );
// Variables for debugging
std::ostringstream os;
if ( Logger::get().level >= Logger::Debug ) {
Debugging( Form( "MaxGen = %d", input_params_->generation.maxgen ) );
}
const unsigned int ndim = function_->dim,
max = pow( mbin_, ndim ),
npoin = input_params_->vegas.npoints;
const double inv_npoin = 1./npoin;
if ( ndim > max_dimensions_ ) {
FatalError( Form( "Number of dimensions to integrate exceed the maximum number, %d", max_dimensions_ ) );
}
nm_ = std::vector<int>( max, 0 );
f_max_ = std::vector<double>( max, 0. );
n_ = std::vector<int>( ndim, 0 );
std::vector<double> x( ndim, 0. );
input_params_->generation.ngen = 0;
// ...
double sum = 0., sum2 = 0., sum2p = 0.;
//--- main loop
for ( unsigned int i=0; i<max; i++ ) {
int jj = i;
for ( unsigned int j=0; j<ndim; j++ ) {
int jjj = jj*inv_mbin_;
n_[j] = jj-jjj*mbin_;
jj = jjj;
}
double fsum = 0., fsum2 = 0.;
for ( unsigned int j=0; j<npoin; j++ ) {
for ( unsigned int k=0; k<ndim; k++ ) {
x[k] = ( uniform() + n_[k] ) * inv_mbin_;
}
const double z = F( x );
f_max_[i] = std::max( f_max_[i], z );
fsum += z;
fsum2 += z*z;
}
const double av = fsum*inv_npoin, av2 = fsum2*inv_npoin, sig2 = av2 - av*av;
sum += av;
sum2 += av2;
sum2p += sig2;
f_max_global_ = std::max( f_max_global_, f_max_[i] );
if ( Logger::get().level >= Logger::DebugInsideLoop ) {
const double sig = sqrt( sig2 );
const double eff = ( f_max_[i] != 0. ) ? f_max_[i]/av : 1.e4;
os.str(""); for ( unsigned int j=0; j<ndim; j++ ) { os << n_[j]; if ( j != ndim-1 ) os << ", "; }
DebuggingInsideLoop( Form( "In iteration #%d:\n\t"
"av = %f\n\t"
"sig = %f\n\t"
"fmax = %f\n\t"
"eff = %f\n\t"
"n = (%s)",
i, av, sig, f_max_[i], eff, os.str().c_str() ) );
}
} // end of main loop
sum = sum/max;
sum2 = sum2/max;
sum2p = sum2p/max;
if ( Logger::get().level >= Logger::Debug ) {
const double sig = sqrt( sum2-sum*sum ), sigp = sqrt( sum2p );
double eff1 = 0.;
for ( unsigned int i=0; i<max; i++ ) eff1 += ( f_max_[i] / ( max*sum ) );
const double eff2 = f_max_global_/sum;
Debugging( Form( "Average function value = sum = %f\n\t"
"Average function value**2 = sum2 = %f\n\t"
"Overall standard deviation = sig = %f\n\t"
"Average standard deviation = sigp = %f\n\t"
"Maximum function value = ffmax = %f\n\t"
"Average inefficiency = eff1 = %f\n\t"
"Overall inefficiency = eff2 = %f\n\t",
sum, sum2, sig, sigp, f_max_global_, eff1, eff2 ) );
}
gen_prepared_ = true;
Information( "Grid prepared! Now launching the production." );
}
}
diff --git a/CepGen/Core/Vegas.h b/CepGen/Core/Vegas.h
index 03bbf1b..eb5e2a6 100644
--- a/CepGen/Core/Vegas.h
+++ b/CepGen/Core/Vegas.h
@@ -1,123 +1,121 @@
#ifndef CepGen_Core_Vegas_h
#define CepGen_Core_Vegas_h
+//#include "CepGen/Parameters.h"
+
#include <gsl/gsl_monte_vegas.h>
#include <gsl/gsl_rng.h>
-#include "CepGen/Parameters.h"
-
#include <vector>
-
-#define fMaxNbins 50
-#define ONE 1.
+#include <memory>
namespace CepGen
{
+ class Parameters;
/**
* Main occurence of the Monte-Carlo integrator @cite PeterLepage1978192 developed by G.P. Lepage in 1978
* \brief Vegas Monte-Carlo integrator instance
*/
class Vegas {
public:
/**
* Book the memory slots and structures for the Vegas integrator
* \note This code is based on the Vegas Monte Carlo integration algorithm developed by P. Lepage, as documented in @cite PeterLepage1978192
* \param[in] dim_ Number of dimensions on which the function will be integrated
* \param[in] f_ Function to be integrated
* \param[inout] inParam_ Run parameters to define the phase space on which this integration is performed (embedded in an Parameters object)
*/
Vegas( const unsigned int dim_, double f_(double*,size_t,void*), Parameters* inParam_ );
/// Class destructor
~Vegas();
/**
* Vegas algorithm to perform the n-dimensional Monte Carlo integration of a given function as described in @cite PeterLepage1978192
* \author Primary author: G.P. Lepage
* \author This C++ implementation: GSL
* \param[out] result_ The cross section as integrated by Vegas for the given phase space restrictions
* \param[out] abserr_ The error associated to the computed cross section
* \return 0 if the integration was performed successfully
*/
int integrate( double& result_,double& abserr_ );
/// Launch the generation of events
void generate();
/**
* Generate one event according to the grid parameters set in Vegas::SetGen
* \brief Generate one single event according to the method defined in the Fortran 77 version of LPAIR
* \return A boolean stating if the generation was successful (in term of the computed weight for the phase space point)
*/
bool generateOneEvent();
/// Dimensional size of the phase space
const unsigned short dimensions() const { return ( !function_ ) ? 0 : function_->dim; }
private:
/**
* Evaluate the function to be integrated at a point @a x_, using the default Parameters object @a fInputParameters
* \param[in] x_ The point at which the function is to be evaluated
* \return Function value at this point @a x_
*/
inline double F( const std::vector<double>& x ) { return F( x, input_params_ ); }
/**
* Evaluate the function to be integrated at a point @a x_, given a set of Parameters @a ip_
* \param[in] x_ The point at which the function is to be evaluated
* \param[in] ip_ A set of parameters to fully define the function
* \return Function value at this point \a x
*/
inline double F( const std::vector<double>& x, Parameters* ip ) {
return function_->f( (double*)&x[0], function_->dim, (void*)ip );
}
/**
* Store the event characterized by its _ndim-dimensional point in the phase
* space to the output file
* \brief Store the event in the output file
* \param[in] x The d-dimensional point in the phase space defining the unique event to store
* \return A boolean stating whether or not the event could be saved
*/
bool storeEvent( const std::vector<double>& x );
/// Start the correction cycle on the grid
/// \param x Point in the phase space considered
/// \param has_correction Correction cycle started?
bool correctionCycle( std::vector<double>& x, bool& has_correction );
/**
* Set all the generation mode variables and align them to the integration grid set while computing the cross-section
* \brief Prepare the class for events generation
*/
void setGen();
double uniform() const { return gsl_rng_uniform( rng_ ); }
- //double uniform() const { return rand()/RAND_MAX; }
/// Maximal number of dimensions handled by this Vegas instance
static constexpr unsigned short max_dimensions_ = 15;
/// Integration grid size parameter
static constexpr unsigned short mbin_ = 3;
static constexpr double inv_mbin_ = 1./mbin_;
/// Selected bin at which the function will be evaluated
int vegas_bin_;
double correc_;
double correc2_;
/// List of parameters to specify the integration range and the physics determining the phase space
Parameters* input_params_;
/// Has the grid been prepared for integration?
bool grid_prepared_;
/// Has the generation been prepared using @a SetGen call? (very time-consuming operation, thus needs to be called once)
bool gen_prepared_;
/// Maximal value of the function at one given point
std::vector<double> f_max_;
double f_max2_;
double f_max_diff_;
double f_max_old_;
/// Maximal value of the function in the considered integration range
double f_max_global_;
std::vector<int> n_;
std::vector<int> nm_;
/// GSL structure storing the function to be integrated by this Vegas instance (along with its parameters)
std::unique_ptr<gsl_monte_function> function_;
gsl_rng* rng_;
/// Number of function calls to be computed for each point
int num_converg_;
/// Number of iterations for the integration
unsigned int num_iter_;
};
}
#endif
diff --git a/CepGen/Core/utils.cpp b/CepGen/Core/utils.cpp
index 7daef12..95a9aae 100644
--- a/CepGen/Core/utils.cpp
+++ b/CepGen/Core/utils.cpp
@@ -1,68 +1,74 @@
-#include "utils.h"
+#include "CepGen/Core/utils.h"
+#include "CepGen/Core/Exception.h"
+
+#include <stdlib.h>
+#include <stdarg.h> // For va_start, etc.
+#include <stdio.h>
+#include <math.h>
void Map( double expo, double xmin, double xmax, double& out, double& dout, const std::string& var_name_ )
{
const double y = xmax/xmin;
out = xmin*pow( y, expo );
dout = out*log( y );
DebuggingInsideLoop( Form( "Mapping variable \"%s\"\n\t"
"min = %f\n\tmax = %f\n\tmax/min = %f\n\t"
"exponent = %f\n\t"
"output = %f\n\td(output) = %f",
var_name_.c_str(), xmin, xmax, y, expo, out, dout ) );
}
void Mapla( double y, double z, int u, double xm, double xp, double& x, double& d )
{
double xmb, xpb, c, yy, zz, alp, alm, am, ap, ax;
xmb = xm-y-z;
xpb = xp-y-z;
c = -4.*y*z;
alp = sqrt( xpb*xpb + c );
alm = sqrt( xmb*xmb + c );
am = xmb+alm;
ap = xpb+alp;
yy = ap/am;
zz = pow( yy, u );
x = y + z + ( am*zz - c / ( am*zz ) ) / 2.;
ax = sqrt( pow( x-y-z, 2 ) + c );
d = ax*log( yy );
}
double BreitWigner( double er, double gamma, double emin, double emax, double x )
{
if ( x==-1. ) x = drand();
if ( gamma<1.e-3*er ) { return er; }
const double a = atan( 2.*( emax-er ) / gamma ),
b = atan( 2.*( emin-er ) / gamma ),
e = er + gamma*tan( x*( a - b ) + b ) / 2.;
if ( e>emax ) { return emax; }
return e;
}
std::string
Form( const std::string fmt, ... )
{
int size = ( (int)fmt.size() ) * 2 + 50; // Use a rubric appropriate for your code
std::string str;
va_list ap;
while ( true ) { // Maximum two passes on a POSIX system...
str.resize( size );
va_start( ap, fmt );
int n = vsnprintf( (char*)str.data(), size, fmt.c_str(), ap );
va_end( ap );
if ( n>-1 and n<size ) { // Everything worked
str.resize( n );
return str;
}
if ( n>-1 ) // Needed size returned
size = n + 1; // For null char
else size *= 2; // Guess at a larger size (OS specific)
}
return str;
}
diff --git a/CepGen/Core/utils.h b/CepGen/Core/utils.h
index 0fcd036..c1820f2 100644
--- a/CepGen/Core/utils.h
+++ b/CepGen/Core/utils.h
@@ -1,55 +1,47 @@
#ifndef CepGen_Core_utils_h
#define CepGen_Core_utils_h
-#include <stdlib.h>
-#include <stdarg.h> // For va_start, etc.
-#include <stdio.h>
-#include <string.h>
-
-#include "CepGen/Core/Exception.h"
-#include "CepGen/Physics/Constants.h"
-
-static std::stringstream error;
+#include <string>
/// Provide a random number generated along a uniform distribution between 0 and 1
//inline double drand() { srand (time(nullptr)); return static_cast<double>(rand())/RAND_MAX; }
#define drand() static_cast<double>( rand()/RAND_MAX )
/// Format a string using a printf style format descriptor.
std::string Form(const std::string fmt, ...);
inline const char* yesno( const bool& test ) { return ( test ) ? "\033[32;1myes\033[0m" : "\033[31;1mno\033[0m"; }
//inline const char* boldify( const char* str ) { const std::string out = std::string( "\033[33;1m" ) + std::string( str ) + std::string( "\033[0m" ); return out.c_str(); }
inline std::string boldify( const std::string& str ) { return Form( "\033[1m%s\033[0m", str.c_str() ); }
inline std::string boldify( const char* str ) { return boldify( std::string( str ) ); }
inline std::string boldify( const double& dbl ) { return boldify( Form("%.2f", dbl ) ); }
inline std::string boldify( const int& i ) { return boldify( Form("% d", i ) ); }
inline std::string boldify( const unsigned int& ui ) { return boldify( Form("%d", ui ) ); }
namespace Colour{
enum TextColour { Gray=30, Red=31, Green=32, Yellow=33, Blue=34, Purple=35 };
}
inline std::string colourise( const std::string& str, const Colour::TextColour& col ) { return Form( "\033[%d%s\033[0m", col, str.c_str() ); }
/**
* Define modified variables of integration to avoid peaks integrations (see @cite Vermaseren1983347 for details)
* Return a set of two modified variables of integration to maintain the stability of the integrant. These two new variables are :
* - \f$y_{out} = x_{min}\left(\frac{x_{max}}{x_{min}}\right)^{exp}\f$ the new variable
* - \f$\mathrm dy_{out} = x_{min}\left(\frac{x_{max}}{x_{min}}\right)^{exp}\log\frac{x_{min}}{x_{max}}\f$, the new variable's differential form
* @brief Redefine the variables of integration in order to avoid the strong peaking of the integrant.
* @param[in] expo Exponant
* @param[in] xmin Minimal value of the variable
* @param[in] xmax Maximal value of the variable
* @param[out] out The new variable definition
* @param[out] dout The differential variant of the new variable definition
* @param[in] var_name The variable name
* @note This method overrides the set of `mapxx` subroutines in ILPAIR, with a slight difference according to the sign of the
* \f$\mathrm dy_{out}\f$ parameter :
* - left unchanged :
* > `mapw2`, `mapxq`, `mapwx`, `maps2`
* - opposite sign :
* > `mapt1`, `mapt2`
*/
void Map( double expo, double xmin, double xmax, double& out, double& dout, const std::string& var_name="" );
void Mapla( double y, double z, int u, double xm, double xp, double& x, double& d );
#endif
diff --git a/CepGen/Event/Event.cpp b/CepGen/Event/Event.cpp
index 9fbc98e..7e853ea 100644
--- a/CepGen/Event/Event.cpp
+++ b/CepGen/Event/Event.cpp
@@ -1,261 +1,262 @@
#include "Event.h"
+#include "CepGen/Core/Exception.h"
#include <algorithm>
namespace CepGen
{
Event::Event() :
num_hadronisation_trials( 0 ),
time_generation( -1. ), time_total( -1. )
{}
Event::~Event()
{}
Event&
Event::operator=( const Event &ev_ )
{
particles_ = ev_.particles_;
time_generation = ev_.time_generation;
time_total = ev_.time_total;
num_hadronisation_trials = ev_.num_hadronisation_trials;
return *this;
}
void
Event::clear()
{
particles_.clear();
time_generation = -1.;
time_total = -1.;
}
void
Event::init()
{
last_particle_ = particles_.end();
}
void
Event::restore()
{
//--- remove all particles after the primordial event block
particles_.erase( last_particle_, particles_.end() );
}
Particles&
Event::getByRole( const Particle::Role& role )
{
//--- retrieve all particles with a given role
return particles_[role];
}
Particle&
Event::getOneByRole( const Particle::Role& role )
{
//--- retrieve the first particle a the given role
Particles& parts_by_role = getByRole( role );
if ( parts_by_role.size() > 1 ) {
FatalError( Form( "More than one particle with role %d: %d particles", (int)role, parts_by_role.size() ) );
}
return *parts_by_role.begin();
}
Particle&
Event::getById( int id )
{
for ( ParticlesMap::iterator out=particles_.begin(); out!=particles_.end(); out++ ) {
for ( Particles::iterator part=out->second.begin(); part!=out->second.end(); part++ ) {
if ( part->id() == id ) return *part;
}
}
throw Exception( __PRETTY_FUNCTION__, Form( "Failed to retrieve the particle with id=%d", id ), FatalError );
}
const Particle&
Event::getConstById( int id ) const
{
for ( ParticlesMap::const_iterator out=particles_.begin(); out!=particles_.end(); out++ ) {
for ( Particles::const_iterator part=out->second.begin(); part!=out->second.end(); part++ ) {
if ( part->id() == id ) return *part;
}
}
throw Exception( __PRETTY_FUNCTION__, Form( "Failed to retrieve the particle with id=%d", id ), FatalError );
}
Particles
Event::getByIds( const ParticlesIds& ids ) const
{
Particles out;
for ( ParticlesIds::const_iterator id=ids.begin(); id!=ids.end(); id++ ) {
out.emplace_back( getConstById( *id ) );
}
return out;
}
Particles
Event::mothers( const Particle& part )
{
return getByIds( part.mothers() );
}
Particles
Event::daughters( const Particle& part )
{
return getByIds( part.daughters() );
}
ParticleRoles
Event::roles() const
{
ParticleRoles out;
ParticlesMap::const_iterator it, end;
for ( it=particles_.begin(), end=particles_.end(); it!=end; it=particles_.upper_bound( it->first ) ) {
out.emplace_back( it->first );
}
return out;
}
void
Event::addParticle( Particle part, bool replace )
{
DebuggingInsideLoop( Form( "Particle with PDGid = %d has role %d", part.pdgId(), part.role() ) );
if ( part.role() <= 0 ) FatalError( Form( "Trying to add a particle with role=%d", (int)part.role() ) );
//--- retrieve the list of particles with the same role
Particles& part_with_same_role = getByRole( part.role() );
//--- specify the id
if ( part_with_same_role.empty() && part.id() < 0 ) part.setId( numParticles() ); // set the id if previously invalid/inexistent
if ( !part_with_same_role.empty() ) {
if ( replace ) part.setId( part_with_same_role[0].id() ); // set the previous id if replacing a particle
else part.setId( numParticles() );
}
//--- add the particle to the collection
if ( replace ) part_with_same_role = Particles( 1, part ); // generate a vector containing only this particle
else part_with_same_role.emplace_back( part );
}
void
Event::addParticle( const Particle::Role& role, bool replace )
{
addParticle( Particle( role ), replace );
}
size_t
Event::numParticles() const
{
size_t out = 0;
for ( ParticlesMap::const_iterator it=particles_.begin(); it!=particles_.end(); it++ ) {
out += it->second.size();
}
return out;
}
const Particles
Event::particles() const
{
Particles out;
for ( ParticlesMap::const_iterator it=particles_.begin(); it!=particles_.end(); it++ ) {
out.insert( out.end(), it->second.begin(), it->second.end() );
}
std::sort( out.begin(), out.end() );
return out;
}
const Particles
Event::stableParticles() const
{
Particles out;
for ( ParticlesMap::const_iterator it=particles_.begin(); it!=particles_.end(); it++ ) {
for ( Particles::const_iterator part=it->second.begin(); part!=it->second.end(); part++ ) {
if ( part->status() == Particle::Undefined || part->status() == Particle::FinalState ) {
out.emplace_back( *part );
}
}
}
std::sort( out.begin(), out.end() );
return out;
}
void
Event::checkKinematics() const
{
try {
const Particles& parts = particles();
// check the kinematics through parentage
for ( Particles::const_iterator p = parts.begin(); p != parts.end(); ++p ) {
ParticlesIds daughters = p->daughters();
if ( daughters.empty() ) continue;
Particle::Momentum ptot;
for ( ParticlesIds::const_iterator daugh = daughters.begin(); daugh != daughters.end(); ++daugh ) {
const Particle& d = getConstById( *daugh );
const ParticlesIds mothers = d.mothers();
if ( mothers.size() > 1 ) {
for ( ParticlesIds::const_iterator moth = mothers.begin(); moth != mothers.end(); ++moth ) {
if ( *moth == p->id() ) continue;
ptot -= getConstById( *moth ).momentum();
}
}
ptot += d.momentum();
}
const double mass_diff = ( ptot-p->momentum() ).mass();
if ( fabs( mass_diff ) > 1.e-10 ) {
dump();
throw Exception( __PRETTY_FUNCTION__, Form( "Error in momentum balance for particle %d: mdiff = %.5e", p->id(), mass_diff ), FatalError );
}
}
} catch ( const Exception& e ) { throw Exception( __PRETTY_FUNCTION__, Form( "Event kinematics check failed:\n\t%s", e.what() ), FatalError ); }
}
void
Event::dump( std::ostream& out, bool stable ) const
{
Particles parts = ( stable ) ? stableParticles() : particles();
std::ostringstream os;
double pxtot = 0., pytot = 0., pztot = 0., etot = 0.;
for ( Particles::const_iterator part_ref=parts.begin(); part_ref!=parts.end(); part_ref++ ) {
const Particle& part = *part_ref;
{ std::ostringstream oss; oss << part.pdgId(); os << Form( "\n %2d\t%+6d%8s", part.id(), part.integerPdgId(), oss.str().c_str() ); }
os << "\t";
if ( part.charge() != 999. ) os << Form( "%6.2f ", part.charge() );
else os << "\t";
{ std::ostringstream oss; oss << part.role(); os << Form( "%8s\t%6d\t", oss.str().c_str(), part.status() ); }
const ParticlesIds mothers = part.mothers();
if ( !mothers.empty() ) {
std::ostringstream oss;
for ( ParticlesIds::const_iterator moth = mothers.begin(); moth != mothers.end(); ++moth ) {
oss << ( ( moth != mothers.begin() ) ? "+" : "" ) << *moth;
}
os << Form( "%6s ", oss.str().c_str() );
}
else os << " ";
const Particle::Momentum mom = part.momentum();
os << Form( "% 9.6e % 9.6e % 9.6e % 9.6e % 9.5e", mom.px(), mom.py(), mom.pz(), part.energy(), part.mass() );
if ( part.status() == Particle::Undefined
|| part.status() == Particle::Undecayed
|| part.status() == Particle::FinalState ) {
const int sign = ( part.status() == Particle::Undefined ) ? -1 : 1;
pxtot += sign*mom.px();
pytot += sign*mom.py();
pztot += sign*mom.pz();
etot += sign*part.energy();
}
}
//--- set a threshold to the computation precision
if ( fabs( pxtot ) < 1.e-10 ) pxtot = 0.;
if ( fabs( pytot ) < 1.e-10 ) pytot = 0.;
if ( fabs( pztot ) < 1.e-10 ) pztot = 0.;
if ( fabs( etot ) < 1.e-10 ) etot = 0.;
//
Information( Form( "Dump of event content:\n"
"Part.\tPDG id\t\tCharge\t Role\tStatus\tMother\t\t\t\t4-Momentum (GeV)\t\tMass (GeV)\n"
"----\t------\t\t------\t ----\t------\t------\t------------------------------------------------------ -----------"
"%s\n"
"---------------------------------------------------------------------------------------------------------------------------\n"
"\t\t\t\t\t\tTotal: % 9.6e % 9.6e % 9.6e % 9.6e", os.str().c_str(), pxtot, pytot, pztot, etot ) );
}
}
diff --git a/CepGen/Event/Event.h b/CepGen/Event/Event.h
index 9096150..9fef87c 100644
--- a/CepGen/Event/Event.h
+++ b/CepGen/Event/Event.h
@@ -1,118 +1,119 @@
#ifndef CepGen_Event_Event_h
#define CepGen_Event_Event_h
#include "Particle.h"
+#include "CepGen/Core/Logger.h"
namespace CepGen
{
/**
* Class containing all the information on the in- and outgoing particles' kinematics
* \brief Kinematic information on the particles in the event
*/
class Event {
public:
Event();
~Event();
/**
* \brief Copies all the relevant quantities from one Event object to another
*/
Event& operator=( const Event& );
/**
* \brief Empties the whole event content
*/
void clear();
/// Initialize an "empty" event collection
void init();
/// Restore the event to its "empty" state
void restore();
/// Dump all the known information on every Particle object contained in this Event container in the output stream
/// \param[out] os Output stream where to dump the information
/// \param[in] stable_ Do we only show the stable particles in this event?
- void dump( std::ostream& os=Logger::get().outputStream, bool stable_=false ) const;
+ void dump( std::ostream& os = Logger::get().outputStream, bool stable_=false ) const;
//----- particles adders
/// Set the information on one particle in the process
/// \param[in] part The Particle object to insert or modify in the event
/// \param[in] replace Do we replace the particle if already present in the event or do we append another particle with the same role ?
void addParticle( Particle part, bool replace=false );
/// \brief Create a new particle in the event, with no kinematic information but the role it has to play in the process
/// \param[in] role The role the particle will play in the process
/// \param[in] replace Do we replace the particle if already present in the event or do we append another particle with the same role ?
void addParticle( const Particle::Role& role, bool replace=false );
//----- particles retrievers
/// Number of particles in the event
size_t numParticles() const;
/// \brief Vector of all particles in the event
const Particles particles() const;
/// \brief Vector of all stable particles in the event
const Particles stableParticles() const;
/**
* Returns the list of Particle objects corresponding to a certain role in the process kinematics
* \brief Gets a list of particles by their role in the event
* \param[in] role The role the particles have to play in the process
* \return A vector of references to the requested Particle objects
*/
Particles& getByRole( const Particle::Role& role );
/**
* Returns the first Particle object in the particles list whose role corresponds to the given argument
* \param[in] role The role the particle has to play in the event
* \return A Particle object corresponding to the first particle with the role
*/
Particle& getOneByRole( const Particle::Role& role );
/**
* Returns the reference to the Particle object corresponding to a unique identifier in the event
* \brief Gets one particle by its unique identifier in the event
* \param[in] id_ The unique identifier to this particle in the event
* \return A reference to the requested Particle object
*/
Particle& getById( int id_ );
/// Get a const Particle object using its unique identifier
/// \param[in] id_ Unique identifier of the particle in the event
/// \return Constant object to be retrieved
const Particle& getConstById( int id_ ) const;
/**
* Returns the references to the Particle objects corresponding to the unique identifiers in the event
* \brief Gets a vector of particles by their unique identifier in the event
* \param[in] ids_ The unique identifiers to the particles to be selected in the event
* \return A vector of references to the requested Particle objects
*/
Particles getByIds( const ParticlesIds& ids_ ) const;
/**
* Returns the list of mother particles of any given Particle object in this event
* \param[in] part The reference to the Particle object from which we want to extract the mother particles
* \return A list of parenting Particle object
*/
//----- general particles information retriever
Particles mothers( const Particle& part );
/// Get a vector containing all the daughters from a particle
/// \param[in] part The particle for which the daughter particles have to be retrieved
/// \return Vector of Particle objects containing all the daughters' kinematic information
Particles daughters( const Particle& part );
/// Get a list of roles for the given event (really process-dependant for the central system)
/// \return Vector of integers corresponding to all the roles the particles can play in the event
ParticleRoles roles() const;
/// Number of trials before the event was "correctly" hadronised
int num_hadronisation_trials;
/// Time needed to generate the event at parton level (in seconds)
float time_generation;
/// Time needed to generate the hadronised (if needed) event (in seconds)
float time_total;
private:
/// Check if the event kinematics is properly defined
void checkKinematics() const;
/// List of particles in the event, mapped to their role in the process
ParticlesMap particles_;
/// Last particle in an "empty" event
ParticlesMap::iterator last_particle_;
};
}
#endif
diff --git a/CepGen/Event/Momentum.cpp b/CepGen/Event/Momentum.cpp
index 03152c0..bcdea69 100644
--- a/CepGen/Event/Momentum.cpp
+++ b/CepGen/Event/Momentum.cpp
@@ -1,286 +1,287 @@
#include "Particle.h"
+#include "CepGen/Core/Exception.h"
namespace CepGen
{
using Momentum = Particle::Momentum;
//----- Particle momentum methods
Momentum::Momentum() :
px_( 0. ), py_( 0. ), pz_( 0. ), p_( 0. ), energy_( 0. )
{}
Momentum::Momentum( double x, double y, double z, double t ) :
px_( x ), py_( y ), pz_( z ), energy_( t )
{
computeP();
}
//--- static constructors
Momentum
Momentum::fromPtEtaPhi( double pt, double eta, double phi, double e )
{
const double px = pt*cos( phi ),
py = pt*sin( phi ),
pz = pt*sinh( eta );
return Momentum( px, py, pz, e );
}
Momentum
Momentum::fromPThetaPhi( double p, double theta, double phi, double e )
{
const double px = p*sin( theta )*cos( phi ),
py = p*sin( theta )*sin( phi ),
pz = p*cos( theta );
return Momentum( px, py, pz, e );
}
Momentum
Momentum::fromPxPyPzE( double px, double py, double pz, double e )
{
return Momentum( px, py, pz, e );
}
//--- arithmetic operators
Momentum&
Momentum::operator+=( const Momentum& mom )
{
px_ += mom.px_;
py_ += mom.py_;
pz_ += mom.pz_;
energy_ += mom.energy_; //FIXME not supposed to be this way!
computeP();
return *this;
}
Momentum&
Momentum::operator-=( const Momentum& mom )
{
px_ -= mom.px_;
py_ -= mom.py_;
pz_ -= mom.pz_;
energy_ -= mom.energy_; //FIXME not supposed to be this way!
computeP();
return *this;
}
double
Momentum::threeProduct( const Momentum& mom ) const
{
DebuggingInsideLoop( Form( " (%f, %f, %f, %f)\n\t* (%f, %f, %f, %f)\n\t= %f",
px_, py_, pz_, energy_,
mom.px_, mom.py_, mom.pz_, mom.energy_,
px_*mom.px_+py_*mom.py_+pz_*mom.pz_ ) );
return px_*mom.px_+py_*mom.py_+pz_*mom.pz_;
}
double
Momentum::fourProduct( const Momentum& mom ) const
{
DebuggingInsideLoop( Form( " (%f, %f, %f, %f)\n\t* (%f, %f, %f, %f)\n\t= %f",
px_, py_, pz_, energy_,
mom.px_, mom.py_, mom.pz_, mom.energy_,
px_*mom.px_+py_*mom.py_+pz_*mom.pz_ ) );
return energy_*mom.energy_-threeProduct(mom);
}
double
Momentum::operator*=( const Momentum& mom )
{
return threeProduct( mom );
}
Momentum&
Momentum::operator*=( double c )
{
px_ *= c;
py_ *= c;
pz_ *= c;
computeP();
return *this;
}
Momentum
operator*( const Momentum& mom, double c )
{
Momentum out = mom;
out *= c;
return out;
}
Momentum
operator*( double c, const Momentum& mom )
{
Momentum out = mom;
out *= c;
return out;
}
Momentum
operator+( const Momentum& mom1, const Momentum& mom2 )
{
Momentum out = mom1;
out += mom2;
return out;
}
Momentum
operator-( const Momentum& mom1, const Momentum& mom2 )
{
Momentum out = mom1;
out -= mom2;
return out;
}
double
operator*( const Momentum& mom1, const Momentum& mom2 )
{
Momentum tmp = mom1;
return tmp.threeProduct( mom2 );
}
//--- various setters
void
Momentum::setP( double px, double py, double pz, double e )
{
setP( px, py, pz );
setEnergy( e );
}
void
Momentum::setP( double px, double py, double pz )
{
px_ = px;
py_ = py;
pz_ = pz;
computeP();
}
void
Momentum::setP( unsigned int i, double p )
{
switch ( i ) {
case 0: px_ = p; break;
case 1: py_ = p; break;
case 2: pz_ = p; break;
case 3: energy_ = p; break;
default: return;
}
computeP();
}
void
Momentum::computeP()
{
p_ = sqrt( px_*px_ + py_*py_ + pz_*pz_ );
}
//--- various getters
double
Momentum::operator[]( const unsigned int i ) const {
switch ( i ) {
case 0: return px_;
case 1: return py_;
case 2: return pz_;
case 3: return energy_;
default: return -1.;
}
}
const std::vector<double>
Momentum::pVector() const
{
return std::vector<double>( { px(), py(), pz(), energy(), mass() } );
}
double
Momentum::mass() const
{
if ( mass2() >= 0. ) return sqrt( mass2() );
return -sqrt( -mass2() );
}
double
Momentum::eta() const
{
const int sign = ( pz()/fabs( pz() ) );
return ( pt() != 0. )
? log( ( p()+fabs( pz() ) )/pt() )*sign
: 9999.*sign;
}
double
Momentum::rapidity() const
{
const int sign = ( pz()/fabs( pz() ) );
return ( energy() >= 0. )
? log( ( energy()+pz() )/( energy()-pz() ) )*0.5
: 999.*sign;
}
//--- boosts/rotations
Momentum&
Momentum::betaGammaBoost( double gamma, double betagamma )
{
if ( gamma == 1. && betagamma == 0. ) return *this; // trivial case
const double pz = pz_, e = energy_;
pz_ = gamma*pz+betagamma*e;
energy_ = gamma*e +betagamma*pz;
computeP();
return *this;
}
Momentum&
Momentum::lorentzBoost( const Momentum& p )
{
const double m = p.mass();
if ( m == p.energy() ) return *this;
const double pf4 = ( ( *this )*p ) / m,
fn = ( pf4+energy() )/( p.energy()+m );
*this -= p*fn;
setEnergy( pf4 );
return *this;
}
Momentum&
Momentum::rotatePhi( double phi, double sign )
{
const double px = px_*cos( phi )+py_*sin( phi )*sign,
py =-px_*sin( phi )+py_*cos( phi )*sign;
px_ = px;
py_ = py;
return *this;
}
Momentum&
Momentum::rotateThetaPhi( double theta, double phi )
{
double rotmtx[3][3], mom[3]; //FIXME check this! cos(phi)->-sin(phi) & sin(phi)->cos(phi) --> phi->phi+pi/2 ?
rotmtx[0][0] = -sin( phi ); rotmtx[0][1] = -cos( theta )*cos( phi ); rotmtx[0][2] = sin( theta )*cos( phi );
rotmtx[1][0] = cos( phi ); rotmtx[1][1] = -cos( theta )*sin( phi ); rotmtx[1][2] = sin( theta )*sin( phi );
rotmtx[2][0] = 0.; rotmtx[2][1] = sin( theta ); rotmtx[2][2] = cos( theta );
for (int i=0; i<3; i++) {
mom[i] = 0.;
for (int j=0; j<3; j++) {
mom[i] += rotmtx[i][j]*operator[]( j );
}
}
setP( mom[0], mom[1], mom[2] );
return *this;
}
//--- printout
std::ostream&
operator<<( std::ostream& os, const Momentum& mom )
{
return os << "(E ; p) = (" << mom.energy_ << " ; " << mom.px_ << ", " << mom.py_ << ", " << mom.pz_ << ")";
}
}
diff --git a/CepGen/Event/Particle.cpp b/CepGen/Event/Particle.cpp
index 20c0c09..27a9241 100644
--- a/CepGen/Event/Particle.cpp
+++ b/CepGen/Event/Particle.cpp
@@ -1,359 +1,361 @@
#include "Particle.h"
+#include "CepGen/Core/Exception.h"
+#include "CepGen/Physics/Constants.h"
namespace CepGen
{
Particle::Particle() :
id_( -1 ), charge_sign_( 1 ),
mass_( -1. ), helicity_( 0. ),
role_( UnknownRole ), status_( Undefined ), pdg_id_( invalidParticle )
{}
Particle::Particle( Role role, ParticleCode pdgId, Status st ) :
id_( -1 ), charge_sign_( 1 ),
mass_( -1. ), helicity_( 0. ),
role_( role ), status_( st ), pdg_id_( pdgId )
{
if ( pdg_id_!=invalidParticle ) {
computeMass();
}
}
Particle::Particle( const Particle& part ) :
id_( part.id_ ), charge_sign_( part.charge_sign_ ),
momentum_( part.momentum_ ), mass_( part.mass_ ), helicity_( part.helicity_ ),
role_( part.role_ ), status_( part.status_ ),
mothers_( part.mothers_ ), daughters_( part.daughters_ ),
pdg_id_( part.pdg_id_ )
{}
bool
Particle::operator<( const Particle& rhs ) const
{
return ( id_ < rhs.id_ );
}
bool
Particle::valid()
{
if ( pdg_id_ == invalidParticle ) return false;
if ( momentum_.p() == 0. && mass() == 0. ) return false;
return true;
}
void
Particle::computeMass( bool off_shell )
{
if ( !off_shell && pdg_id_ != invalidParticle ) { // retrieve the mass from the on-shell particle's properties
mass_ = massFromPDGId( pdg_id_ );
}
else if ( momentum_.energy() >= 0. ) {
mass_ = sqrt( energy2() - momentum_.p2() );
}
//--- finish by setting the energy accordingly
if ( momentum_.energy() < 0. ) { // invalid energy
momentum_.setEnergy( sqrt( momentum_.p2() + mass2() ) );
}
}
void
Particle::setMass( double m )
{
if ( m >= 0. ) mass_ = m;
else computeMass();
}
void
Particle::addMother( Particle& part )
{
mothers_.insert( part.id() );
DebuggingInsideLoop( Form( "Particle %2d (pdgId=%4d) is the new mother of %2d (pdgId=%4d)",
part.id(), part.pdgId(), id_, pdg_id_ ) );
part.addDaughter( *this );
}
void
Particle::addDaughter( Particle& part )
{
std::pair<ParticlesIds::iterator,bool> ret = daughters_.insert( part.id() );
if ( Logger::get().level >= Logger::DebugInsideLoop ) {
std::ostringstream os;
for ( ParticlesIds::const_iterator it=daughters_.begin(); it!=daughters_.end(); it++) {
os << Form( "\n\t * id=%d", *it );
}
DebuggingInsideLoop( Form( "Particle %2d (pdgId=%4d) has now %2d daughter(s):"
"%s", role_, pdg_id_, numDaughters(), os.str().c_str() ) );
}
if ( ret.second ) {
DebuggingInsideLoop( Form( "Particle %2d (pdgId=%4d) is a new daughter of %2d (pdgId=%4d)",
part.role(), part.pdgId(), role_, pdg_id_ ) );
if ( part.mothers().find( id_ ) == part.mothers().end() ) part.addMother( *this );
}
}
void
Particle::setMomentum( const Momentum& mom, bool offshell )
{
momentum_ = mom;
if ( !offshell && mom.mass() > 0. ) mass_ = momentum_.mass();
else computeMass();
}
void
Particle::setMomentum( double px, double py, double pz )
{
momentum_.setP( px, py, pz );
setEnergy();
}
void
Particle::setMomentum( double px, double py, double pz, double e )
{
setMomentum( px, py, pz );
if ( fabs( e-momentum_.energy() )>1.e-6 ) { // more than 1 eV difference
InError( Form( "Energy difference: %.5e", e-momentum_.energy() ) );
return;
}
}
void
Particle::setEnergy( double e )
{
if ( e < 0. && mass_ >= 0. ) e = sqrt( mass2()+momentum_.p2() );
momentum_.setEnergy( e );
}
void
Particle::setPdgId( const ParticleCode& pdg, short ch )
{
pdg_id_ = pdg;
charge_sign_ = ch;
}
int
Particle::integerPdgId() const
{
const float ch = chargeFromPDGId( pdg_id_ );
if ( ch == 0 ) return static_cast<int>( pdg_id_ );
return static_cast<int>( pdg_id_ ) * charge_sign_ * ( ch/fabs( ch ) );
}
void
Particle::dump() const
{
std::ostringstream osm, osd, os;
if ( !primary() ) {
osm << ": mother(s): ";
for ( ParticlesIds::const_iterator m=mothers_.begin(); m!=mothers_.end(); m++ ) {
if ( m!=mothers_.begin() ) osm << ", ";
osm << ( *m );
}
}
const ParticlesIds daugh = daughters();
if ( daugh.size()!=0 ) {
osd << ": id = ";
for ( ParticlesIds::const_iterator it=daugh.begin(); it!=daugh.end(); it++ ) {
if ( it!=daugh.begin() ) osd << ", ";
osd << ( *it );
}
}
os << " (" << pdg_id_ << ")";
if ( os.str() == " ()" ) os.str("");
Information( Form(
"Dumping a particle with id=%3d, role=%3d, status=% 3d\n\t"
"PDG Id:%4d%s, mass = %5.4f GeV\n\t"
"(E,P) = (%4.2f, %4.2f, %4.2f, %4.2f) GeV\t"
"(|P| = p = %4.2f GeV)\n\t"
" Pt = %5.4f GeV, eta = %4.3f, phi = % 4.3f\n\t"
"Primary? %s%s\n\t"
"%d daughter(s)%s",
id_, role_, status_, pdg_id_, os.str().c_str(),
mass(), energy(), momentum_.px(), momentum_.py(), momentum_.pz(),
momentum_.p(), momentum_.pt(), momentum_.eta(), momentum_.phi(),
yesno( primary() ), osm.str().c_str(), numDaughters(), osd.str().c_str() )
);
}
Particle&
Particle::lorentzBoost( double m, const Particle::Momentum& mom )
{
double pf4, fn;
if ( mom.energy() != m ) {
pf4 = 0.;
for ( unsigned int i=0; i<4; i++ ) {
pf4 += momentum_[i]*mom[i];
}
pf4 /= m;
fn = ( pf4+energy() )/( momentum_.energy()+m );
for ( unsigned int i=0; i<3; i++ ) {
momentum_.setP( i, momentum_[i]+fn*mom[i] );
}
}
return *this;
}
std::vector<double>
Particle::lorentzBoost( const Particle::Momentum& mom )
{
double p2, gamma, bp, gamma2;
std::vector<double> out( 3, 0. );
p2 = mom.p2();
gamma = 1./sqrt( 1.-p2 );
bp = 0.;
for ( unsigned int i=0; i<3; i++ ) bp+= mom[i]*momentum_[i];
if ( p2>0. ) gamma2 = (gamma-1.)/p2;
else gamma2 = 0.;
for ( unsigned int i=0; i<3; i++ ) {
out[i] = momentum_[i] + gamma2*bp*mom[i]+gamma*mom[i]*energy();
}
return out;
}
double
Particle::massFromPDGId( const Particle::ParticleCode& pdg )
{
switch ( pdg ) {
case dQuark: return 0.33; // mass from PYTHIA6.4
case uQuark: return 0.33; // mass from PYTHIA6.4
case Electron: return 0.510998928e-3;
case ElectronNeutrino: return 0.;
case Muon: return 0.1056583715;
case MuonNeutrino: return 0.;
case Tau: return 1.77682;
case TauNeutrino: return 0.;
case Gluon: return 0.;
case Z: return 91.1876;
case W: return 80.385;
case Photon: return 0.;
case PiPlus: return 0.13957018;
case PiZero: return Constants::mpi;
case JPsi: return 20.; //FIXME FIXME FIXME
case ud0Diquark: return 0.57933;
case ud1Diquark: return 0.77133;
case uu1Diquark: return 0.77133;
case Proton: return Constants::mp;
case Neutron: return 0.939565346;
case Upsilon1S: return 9.46030;
case Upsilon2S: return 10.02326;
case Upsilon3S: return 10.3552;
case Rho770_0: return 0.77526;
case Rho1450_0: return 1.465;
case Rho1700_0: return 1.720;
case h1380_1: return 1.38619;
case Eta: return 0.547862;
case invalidParticle:
default: return -1.;
}
}
double
Particle::chargeFromPDGId( const Particle::ParticleCode& pdg )
{
switch ( pdg ) {
case Proton: return +1.;
case dQuark: return -1./3;
case uQuark: return +2./3;
case Electron: case Muon: case Tau: return -1.;
case ElectronNeutrino: case MuonNeutrino: case TauNeutrino: return 0.;
case Gluon: case Z: case Photon: return 0.;
case W: return +1.;
case PiPlus: return +1.;
case PiZero: return 0.;
case Neutron: return 0.;
case Eta: return 0.;
default: return 0.;
}
}
double
Particle::widthFromPDGId( const Particle::ParticleCode& pdg )
{
switch ( pdg ) {
case JPsi: return 5.; //FIXME
case Z: return 2.4952;
case W: return 2.085;
case Upsilon1S: return 54.02e-6;
case Upsilon2S: return 31.98e-6;
case Upsilon3S: return 20.32e-6;
case Rho770_0: return 0.150; // PDG
case Rho1450_0: return 0.400; // PDG
case Rho1700_0: return 0.250; // PDG
default: return -1.;
}
}
double
Particle::etaToY( double eta_, double m_, double pt_ )
{
const double mt = m_*m_ + pt_*pt_;
return asinh( sqrt( ( ( ( mt*mt-m_*m_ )*cosh( 2.*eta_ ) + m_*m_ )/ mt*mt - 1. ) / 2. ) );
}
std::ostream&
operator<<( std::ostream& os, const Particle::ParticleCode& pc )
{
switch ( pc ) {
case Particle::dQuark: return os << "d quark";
case Particle::uQuark: return os << "u quark";
case Particle::Electron: return os << "electron";
case Particle::ElectronNeutrino: return os << "nu_e";
case Particle::Muon: return os << "muon";
case Particle::MuonNeutrino: return os << "nu_mu";
case Particle::Tau: return os << "tau";
case Particle::TauNeutrino: return os << "nu_tau";
case Particle::Gluon: return os << "gluon";
case Particle::Photon: return os << "photon";
case Particle::Z: return os << "Z";
case Particle::W: return os << "W+-";
case Particle::PiPlus: return os << "pi+";
case Particle::PiZero: return os << "pi0";
case Particle::Rho770_0: return os << "rho(770)0";
case Particle::Rho1450_0: return os << "rho(1450)0";
case Particle::Rho1700_0: return os << "rho(1700)0";
case Particle::h1380_1: return os << "h(1380)1";
case Particle::Eta: return os << "eta meson";
case Particle::Omega782: return os << "omega(782)";
case Particle::JPsi: return os << "J/Psi";
case Particle::Phi1680: return os << "phi(1680)";
case Particle::Upsilon1S: return os << "Upsilon(1S)";
case Particle::Upsilon2S: return os << "Upsilon(2S)";
case Particle::Upsilon3S: return os << "Upsilon(3S)";;
case Particle::ud0Diquark: return os << "(ud)0 di-quark";
case Particle::ud1Diquark: return os << "(ud)1 di-quark";
case Particle::uu1Diquark: return os << "(uu)1 di-quark";
case Particle::Proton: return os << "proton";
case Particle::Neutron: return os << "neutron";
case Particle::Pomeron: return os << "pomeron";
case Particle::Reggeon: return os << "reggeon";
case Particle::invalidParticle: return os << "[...]";
}
return os;
}
std::ostream&
operator<<( std::ostream& os, const Particle::Role& rl )
{
switch ( rl ) {
case Particle::UnknownRole: return os << "unknown";
case Particle::IncomingBeam1: return os << "in.b.1";
case Particle::IncomingBeam2: return os << "in.b.2";
case Particle::OutgoingBeam1: return os << "out.b.1";
case Particle::OutgoingBeam2: return os << "out.b.2";
case Particle::Parton1: return os << "parton1";
case Particle::Parton2: return os << "parton2";
case Particle::Parton3: return os << "parton3";
case Particle::Intermediate: return os << "partons";
case Particle::CentralSystem: return os << "central";
}
return os;
}
}
diff --git a/CepGen/Event/Particle.h b/CepGen/Event/Particle.h
index db1901d..133466b 100644
--- a/CepGen/Event/Particle.h
+++ b/CepGen/Event/Particle.h
@@ -1,415 +1,416 @@
#ifndef CepGen_Event_Particle_h
#define CepGen_Event_Particle_h
#include <set>
#include <map>
#include <vector>
#include "CepGen/Core/utils.h"
+#include <cmath>
namespace CepGen
{
/// A set of integer-type particle identifiers
typedef std::set<int> ParticlesIds;
/// Kinematic information for one particle
class Particle {
public:
/** Unique identifier for a particle type. From \cite Beringer:1900zz :
* `The Monte Carlo particle numbering scheme [...] is intended to facilitate interfacing between event generators, detector simulators, and analysis packages used in particle physics.`
* \brief PDG ids of all known particles
*/
enum ParticleCode {
invalidParticle = 0,
//--- fundamental particles
dQuark = 1, uQuark = 2,
Electron = 11, ElectronNeutrino = 12,
Muon = 13, MuonNeutrino = 14,
Tau = 15, TauNeutrino = 16,
Gluon = 21, Photon = 22, Z = 23, W = 24,
//--- composite particles
PiPlus = 211, PiZero = 111,
Rho770_0 = 113, Rho1450_0 = 100113, Rho1700_0 = 30113,
Eta = 221, Omega782 = 223,
h1380_1 = 10333,
JPsi= 443,
Phi1680 = 100333,
Upsilon1S = 553, Upsilon2S = 100553, Upsilon3S = 200553,
ud0Diquark = 2101, ud1Diquark = 2103, uu1Diquark = 2203,
Proton = 2212, Neutron = 2112,
Pomeron = 990, Reggeon = 110
};
/// Internal status code for a particle
enum Status {
PrimordialIncoming = -9,
Undecayed = -3,
Propagator = -2,
Incoming = -1,
Undefined = 0,
FinalState = 1,
Resonance = 2,
DebugResonance = 3,
PythiaHIncoming = 21
};
/// Role of the particle in the process
enum Role {
UnknownRole = -1,
IncomingBeam1 = 1, IncomingBeam2 = 2,
OutgoingBeam1 = 3, OutgoingBeam2 = 5,
CentralSystem = 6,
Intermediate = 4,
Parton1 = 41, Parton2 = 42, Parton3 = 43
};
/**
* Container for a particle's 4-momentum, along with useful methods to ease the development of any matrix element level generator
* \brief 4-momentum for a particle
* \date Dec 2015
* \author Laurent Forthomme <laurent.forthomme@cern.ch>
*/
class Momentum {
public:
/// Build a 4-momentum at rest with an invalid energy (no mass information known)
Momentum();
/// Build a 4-momentum using its 3-momentum coordinates and its energy
Momentum( double x_, double y_, double z_, double t_=-1. );
inline ~Momentum() {}
// --- static definitions
/// Build a 3-momentum from its three pseudo-cylindric coordinates
static Momentum fromPtEtaPhi( double pt, double eta, double phi, double e=-1. );
/// Build a 4-momentum from its scalar momentum, and its polar and azimuthal angles
static Momentum fromPThetaPhi( double p, double theta, double phi, double e=-1. );
/// Build a 4-momentum from its four momentum and energy coordinates
static Momentum fromPxPyPzE( double px, double py, double pz, double e );
// --- vector and scalar operators
/// Scalar product of the 3-momentum with another 3-momentum
double threeProduct( const Momentum& ) const;
/// Scalar product of the 4-momentum with another 4-momentum
double fourProduct( const Momentum& ) const;
/// Add a 4-momentum through a 4-vector sum
Momentum& operator+=( const Momentum& );
/// Subtract a 4-momentum through a 4-vector sum
Momentum& operator-=( const Momentum& );
/// Scalar product of the 3-momentum with another 3-momentum
double operator*=( const Momentum& );
/// Multiply all 4-momentum coordinates by a scalar
Momentum& operator*=( double c );
/// Human-readable format for a particle's momentum
friend std::ostream& operator<<( std::ostream& os, const Particle::Momentum& mom );
Momentum& betaGammaBoost( double gamma, double betagamma );
/// Forward Lorentz boost
Momentum& lorentzBoost( const Particle::Momentum& p );
// --- setters and getters
/// Set all the components of the 4-momentum (in GeV)
void setP( double px, double py, double pz, double e );
/// Set all the components of the 3-momentum (in GeV)
void setP( double px, double py, double pz );
/// Set an individual component of the 4-momentum (in GeV)
void setP( unsigned int i, double p );
/// Set the energy (in GeV)
inline void setEnergy( double e ) { energy_ = e; }
/// Compute the energy from the mass
inline void setMass( double m ) { energy_ = sqrt( p2()+m*m ); }
/// Compute the energy from the mass
inline void setMass2( double m2 ) { energy_ = sqrt( p2()+m2 ); }
/// Get one component of the 4-momentum (in GeV)
double operator[]( const unsigned int i ) const;
/// Momentum along the \f$x\f$-axis (in GeV)
inline double px() const { return px_; }
/// Momentum along the \f$y\f$-axis (in GeV)
inline double py() const { return py_; }
/// Longitudinal momentum (in GeV)
inline double pz() const { return pz_; }
/// Transverse momentum (in GeV)
inline double pt() const { return sqrt( pt2() ); }
/// Squared transverse momentum (in GeV\f$^\textrm{2}\f$)
inline double pt2() const { return ( px()*px()+py()*py() ); }
/// 4-vector of double precision floats (in GeV)
const std::vector<double> pVector() const;
/// 3-momentum norm (in GeV)
inline double p() const { return p_; }
/// Squared 3-momentum norm (in GeV\f$^\textrm{2}\f$)
inline double p2() const { return p_*p_; }
/// Energy (in GeV)
inline double energy() const { return energy_; }
/// Squared energy (in GeV^2)
inline double energy2() const { return energy_*energy_; }
/// Squared mass (in GeV^2) as computed from its energy and momentum
inline double mass2() const { return energy2()-p2(); }
/// Mass (in GeV) as computed from its energy and momentum
/// \note Returns \f$-\sqrt{|E^2-\mathbf{p}^2|}<0\f$ if \f$\mathbf{p}^2>E^2\f$
double mass() const;
/// Polar angle (angle with respect to the longitudinal direction)
inline double theta() const { return atan2( pt(), pz() ); }
/// Azimutal angle (angle in the transverse plane)
inline double phi() const { return atan2( py(), px() ); }
/// Pseudo-rapidity
double eta() const;
/// Rapidity
double rapidity() const;
/// Rotate the transverse components by an angle phi (and reflect the y coordinate)
Momentum& rotatePhi( double phi, double sign );
/// Rotate the particle's momentum by a polar/azimuthal angle
Momentum& rotateThetaPhi( double theta_, double phi_ );
/// Apply a \f$ z\rightarrow -z\f$ transformation
inline Momentum& mirrorZ() { pz_ = -pz_; return *this; }
private:
/// Compute the 3-momentum's norm
void computeP();
/// Momentum along the \f$x\f$-axis
double px_;
/// Momentum along the \f$y\f$-axis
double py_;
/// Momentum along the \f$z\f$-axis
double pz_;
/// 3-momentum's norm (in GeV/c)
double p_;
/// Energy (in GeV)
double energy_;
};
/// Human-readable format for a particle's PDG code
friend std::ostream& operator<<( std::ostream& os, const Particle::ParticleCode& pc );
/// Human-readable format for a particle's role in the event
friend std::ostream& operator<<( std::ostream& os, const Particle::Role& rl );
/// Compute the 4-vector sum of two 4-momenta
friend Particle::Momentum operator+( const Particle::Momentum& mom1, const Particle::Momentum& mom2 );
/// Compute the 4-vector difference of two 4-momenta
friend Particle::Momentum operator-( const Particle::Momentum& mom1, const Particle::Momentum& mom2 );
/// Scalar product of two 3-momenta
friend double operator*( const Particle::Momentum& mom1, const Particle::Momentum& mom2 );
/// Multiply all components of a 4-momentum by a scalar
friend Particle::Momentum operator*( const Particle::Momentum& mom, double c );
/// Multiply all components of a 4-momentum by a scalar
friend Particle::Momentum operator*( double c, const Particle::Momentum& mom );
//----- static getters
/**
* \brief Gets the mass of a particle
* \param pdgId ParticleCode (PDG ID)
* \return Mass of the particle in \f$\textrm{GeV}/c^2\f$
*/
static double massFromPDGId( const Particle::ParticleCode& pdgId );
/**
* \brief Gets the electric charge of a particle
* \param pdgId ParticleCode (PDG ID)
* \return Charge of the particle in \f$e\f$
*/
static double chargeFromPDGId( const Particle::ParticleCode& pdgId );
/**
* \brief Total decay width of one unstable particle
* \param[in] pdgId ParticleCode (PDG ID)
* \return Decay width in GeV
*/
static double widthFromPDGId( const Particle::ParticleCode& pdgId );
/// Convert a polar angle to a pseudo-rapidity
static inline double thetaToEta( double theta_ ) { return -log( tan( theta_/180.*M_PI/2. ) ); }
/// Convert a pseudo-rapidity to a polar angle
static inline double etaToTheta( double eta_ ) { return 2.*atan( exp( -eta_ ) )*180. / M_PI; }
/// Convert a pseudo-rapidity to a rapidity
static double etaToY( double eta_, double m_, double pt_ );
Particle();
/// Build using the role of the particle in the process and its PDG id
/// \param[in] pdgId ParticleCode (PDG ID)
/// \param[in] role Role of the particle in the process
/// \param[in] st Current status
Particle( Role role, ParticleCode pdgId=invalidParticle, Status st=Undefined );
/// Copy constructor
Particle( const Particle& );
inline ~Particle() {}
/// Comparison operator (from unique identifier)
bool operator<( const Particle& rhs ) const;
/// Comparison operator (from their reference's unique identifier)
//bool operator<( Particle *rhs ) const { return ( id < rhs->id ); }
Particle& lorentzBoost( double m_, const Momentum& mom_ );
/// Lorentz boost (shamelessly stolen from ROOT)
std::vector<double> lorentzBoost( const Momentum& mom_ );
// --- general particle properties
/// Unique identifier (in a Event object context)
int id() const { return id_; }
//void setId( int id ) { id_ = id; }
/// Set the particle unique identifier in an event
void setId( int id ) { id_ = id; }
/// Electric charge (given as a float number, for the quarks and bound states)
float charge() const { return charge_sign_ * chargeFromPDGId( pdg_id_ ); }
/// Set the electric charge sign (+-1 for charged or 0 for neutral particles)
void setChargeSign( int sign ) { charge_sign_ = sign; }
/// Role in the considered process
Role role() const { return role_; }
/// Set the particle role in the process
void setRole( const Role& role ) { role_ = role; }
/**
* Codes 1-10 correspond to currently existing partons/particles, and larger codes contain partons/particles which no longer exist, or other kinds of event information
* \brief Particle status
*/
Status status() const { return status_; }
/// Set the particle decay/stability status
void setStatus( const Status& status ) { status_ = status; }
/// Set the PDG identifier (along with the particle's electric charge)
/// \param[in] pdg ParticleCode (PDG ID)
/// \param[in] ch Electric charge (0, 1, or -1)
void setPdgId( const ParticleCode& pdg, short ch=0 );
/// Retrieve the objectified PDG identifier
inline ParticleCode pdgId() const { return pdg_id_; }
/// Retrieve the integer value of the PDG identifier
int integerPdgId() const;
/// Particle's helicity
float helicity() const { return helicity_; }
/// Set the helicity of the particle
void setHelicity( float heli ) { helicity_ = heli; }
/**
* Gets the particle's mass in \f$\textrm{GeV}/c^{2}\f$.
* \brief Gets the particle's mass
* \return The particle's mass
*/
inline double mass() const { return mass_; };
/**
* Set the mass of the particle in \f$\textrm{GeV}/c^{2}\f$ while ensuring that the kinematics is properly set (the mass is set according to the energy and the momentum in priority)
* \brief Compute the particle's mass in \f$\textrm{GeV}/c^{2}\f$
*/
void computeMass( bool off_shell=false );
/**
* Set the mass of the particle in \f$\textrm{GeV}/c^{2}\f$ according to a value given as an argument. This method ensures that the kinematics is properly set (the mass is set according to the energy and the momentum in priority)
* \param m The mass in \f$\textrm{GeV}/c^{2}\f$ to set
* \brief Set the particle's mass in \f$\textrm{GeV}/c^{2}\f$
*/
void setMass( double m=-1. );
/// Get the particle's squared mass (in \f$\textrm{GeV}^\textrm{2}\f$)
inline double mass2() const { return mass_*mass_; };
/// Retrieve the momentum object associated with this particle
inline Momentum& momentum() { return momentum_; }
/// Retrieve the momentum object associated with this particle
inline Momentum momentum() const { return momentum_; }
/// Associate a momentum object to this particle
void setMomentum( const Momentum& mom, bool offshell=false );
/**
* \brief Set the 3-momentum associated to the particle
* \param[in] px Momentum along the \f$x\f$-axis, in \f$\textrm{GeV}/c\f$
* \param[in] py Momentum along the \f$y\f$-axis, in \f$\textrm{GeV}/c\f$
* \param[in] pz Momentum along the \f$z\f$-axis, in \f$\textrm{GeV}/c\f$
*/
void setMomentum( double px, double py, double pz );
/**
* \brief Set the 4-momentum associated to the particle
* \param[in] px Momentum along the \f$x\f$-axis, in \f$\textrm{GeV}/c\f$
* \param[in] py Momentum along the \f$y\f$-axis, in \f$\textrm{GeV}/c\f$
* \param[in] pz Momentum along the \f$z\f$-axis, in \f$\textrm{GeV}/c\f$
* \param[in] e Energy, in GeV
*/
void setMomentum( double px, double py, double pz, double e );
/**
* \brief Set the 4-momentum associated to the particle
* \param[in] p 4-momentum
*/
inline void setMomentum( double p[4] ) { setMomentum( p[0], p[1], p[2], p[3] ); }
/**
* \brief Set the particle's energy
* \param[in] e Energy, in GeV
*/
void setEnergy( double e=-1. );
/// Get the particle's energy (in GeV)
inline double energy() const {
return ( momentum_.energy() < 0. ) ? sqrt( mass2()+momentum_.p2() ) : momentum_.energy();
};
/// Get the particle's squared energy (in \f$\textrm{GeV}^\textrm{2}\f$)
inline double energy2() const { return energy()*energy(); };
/// Is this particle a valid particle which can be used for kinematic computations?
bool valid();
// --- particle relations
/// Is this particle a primary particle?
inline bool primary() const { return mothers_.empty(); }
/**
* \brief Set the mother particle
* \param[in] part A Particle object containing all the information on the mother particle
*/
void addMother( Particle& part );
/**
* \brief Gets the unique identifier to the mother particle from which this particle arises
* \return An integer representing the unique identifier to the mother of this particle in the event
*/
inline ParticlesIds mothers() const { return mothers_; }
/**
* \brief Add a decay product
* \param[in] part The Particle object in which this particle will desintegrate or convert
* \return A boolean stating if the particle has been added to the daughters list or if it was already present before
*/
void addDaughter( Particle& part );
/// Gets the number of daughter particles
inline unsigned int numDaughters() const { return daughters_.size(); };
/**
* \brief Get an identifiers list all daughter particles
* \return An integer vector containing all the daughters' unique identifier in the event
*/
inline ParticlesIds daughters() const { return daughters_; }
// --- global particle information extraction
/// Dump all the information on this particle into the standard output stream
void dump() const;
private:
/// Unique identifier in an event
int id_;
/// Electric charge (+-1 or 0)
short charge_sign_;
/// Momentum properties handler
Momentum momentum_;
/// Mass in \f$\textrm{GeV}/c^2\f$
double mass_;
/// Helicity
float helicity_;
/// Role in the process
Role role_;
/// Decay/stability status
Status status_;
/// List of mother particles
ParticlesIds mothers_;
/// List of daughter particles
ParticlesIds daughters_;
/// PDG id
ParticleCode pdg_id_;
};
/// Compute the centre of mass energy of two particles (incoming or outgoing states)
inline static double CMEnergy( const Particle& p1, const Particle& p2 ) {
if ( p1.mass()*p2.mass() < 0. ) return 0.;
if ( p1.energy()*p2.energy() < 0. ) return 0.;
return sqrt( p1.mass2()+p2.mass2() + 2.*p1.energy()*p2.energy() - 2.*( p1.momentum()*p2.momentum() ) );
}
/// Compute the centre of mass energy of two particles (incoming or outgoing states)
inline static double CMEnergy( const Particle::Momentum& m1, const Particle::Momentum& m2 ) {
if ( m1.mass()*m2.mass() < 0. ) return 0.;
if ( m1.energy()*m2.energy() < 0. ) return 0.;
return sqrt( m1.mass2()+m2.mass2() + 2.*m1.energy()*m2.energy() - 2.*( m1*m2 ) );
}
//bool operator<( const Particle& a, const Particle& b ) { return a.id<b.id; }
// --- particle containers
/// List of Particle objects
typedef std::vector<Particle> Particles;
/// List of particles' roles
typedef std::vector<Particle::Role> ParticleRoles;
/// Map between a particle's role and its associated Particle object
typedef std::map<Particle::Role,Particles> ParticlesMap;
}
#endif
diff --git a/CepGen/Generator.h b/CepGen/Generator.h
index 5a5de21..5f7672f 100644
--- a/CepGen/Generator.h
+++ b/CepGen/Generator.h
@@ -1,131 +1,123 @@
#ifndef CepGen_Generator_h
#define CepGen_Generator_h
#include <sstream>
#include <memory>
-#include "CepGen/Core/Vegas.h"
-
////////////////////////////////////////////////////////////////////////////////
/**
* \image latex cepgen_logo.pdf
* \mainpage Foreword
* This Monte Carlo generator was developed as a modern version of the LPAIR code introduced
* in the early 1990s by J. Vermaseren *et al*\cite Vermaseren1983347. This latter allows to
* compute the cross-section and to generate events for the \f$\gamma\gamma\to\ell^{+}\ell^{-}\f$
* process in the scope of high energy physics.
*
* Soon after the integration of its matrix element, it was extended as a tool to compute and
* generate events for any generic 2\f$\rightarrow\f$ 3 central exclusive process.
* To do so, the main operation performed here is the integration of the matrix element (given as a
* subset of a GenericProcess object) by the GSL implementation of the *Vegas* algorithm, a
* numerical technique for importance sampling integration developed in 1972 by G. P. Lepage\cite PeterLepage1978192.
*
*/
////////////////////////////////////////////////////////////////////////////////
/// Common namespace for this Monte Carlo generator
namespace CepGen
{
/**
* Function to be integrated. It returns the value of the weight for one point
* of the full phase space (or "event"). This weights includes the matrix element
* of the process considered, along with all the kinematic factors, and the cut
* restrictions imposed on this phase space. \f$x\f$ is therefore an array of random
* numbers defined inside its boundaries (as normalised so that \f$\forall i<\mathrm{ndim}\f$,
* \f$0<x_i<1\f$.
*/
double f( double*, size_t, void* );
+ class Event;
+ class Vegas;
+ class Parameters;
+
////////////////////////////////////////////////////////////////////////////////
/**
* This object represents the core of this Monte Carlo generator, with its
* capability to generate the events (using the embedded Vegas object) and to
* study the phase space in term of the variation of resulting cross section
* while scanning the various parameters (point \f$\textbf{x}\f$ in the
* multi-dimensional phase space).
*
* The phase space is constrained using the Parameters object given as an
* argument to the constructor, and the differential cross-sections for each
* value of the array \f$\textbf{x}\f$ are computed in the \a f-function defined
* outside (but populated inside) this object.
*
* This f-function embeds a GenericProcess-inherited object which defines all the
* methods to compute this differential cross-section as well as the in- and outgoing
* kinematics associated to each particle.
*
* \author Laurent Forthomme <laurent.forthomme@cern.ch>
* \date Feb 2013
* \brief Core of the Monte-Carlo generator
*
*/
class Generator {
public:
/// Core of the Monte Carlo integrator and events generator
Generator();
/// Core of the Monte Carlo integrator and events generator
/// \param[in] ip List of input parameters defining the phase space on which to perform the integration
Generator( Parameters *ip );
~Generator();
/// Dump this program's header into the standard output stream
void printHeader();
/// Feed the generator with a Parameters object
void setParameters( Parameters& ip );
/// Remove all references to a previous generation/run
void clearRun();
/**
* Compute the cross section for the run parameters defined by this object.
* This returns the cross section as well as the absolute error computed along.
* \brief Compute the cross-section for the given process
* \param[out] xsec The computed cross-section, in pb
* \param[out] err The absolute integration error on the computed cross-section, in pb
*/
void computeXsection( double& xsec, double& err );
/// Last cross section computed by the generator
double crossSection() const { return cross_section_; }
/// Last error on the cross section computed by the generator
double crossSectionError() const { return cross_section_error_; }
/**
* Generate one single event given the phase space computed by Vegas in the integration step
* \return A pointer to the Event object generated in this run
*/
Event* generateOneEvent();
/// Number of dimensions on which the integration is performed
- inline size_t numDimensions() const {
- if ( !parameters->process() ) return 0;
- return parameters->process()->numDimensions( parameters->kinematics.mode );
- }
+ size_t numDimensions() const;
/// Compute one single point from the total phase space
- /// \param[in] x_ the n-dimensional point to compute
+ /// \param[in] x the n-dimensional point to compute
/// \return the function value for the given point
- inline double computePoint( double* x_ ) {
- prepareFunction();
- double res = f( x_, numDimensions(), (void*)parameters.get() );
- std::ostringstream os;
- for ( unsigned int i=0; i<numDimensions(); i++ ) { os << x_[i] << " "; }
- Debugging( Form( "Result for x[%zu] = ( %s):\n\t%10.6f", numDimensions(), os.str().c_str(), res ) );
- return res;
- }
+ double computePoint( double* x );
/// Physical Parameters used in the events generation and cross-section computation
std::unique_ptr<Parameters> parameters;
/// Last event generated in this run
std::shared_ptr<Event> last_event;
private:
/// Prepare the function before its integration (add particles/compute kinematics/...)
void prepareFunction();
/// Vegas instance which will integrate the function
std::unique_ptr<Vegas> vegas_;
/// Cross section value computed at the last integration
double cross_section_;
/// Error on the cross section as computed in the last integration
double cross_section_error_;
/// Has a first integration beed already performed?
bool has_cross_section_;
};
}
#endif
diff --git a/CepGen/Physics/FormFactors.cpp b/CepGen/Physics/FormFactors.cpp
index e8c5909..c92f0bd 100644
--- a/CepGen/Physics/FormFactors.cpp
+++ b/CepGen/Physics/FormFactors.cpp
@@ -1,76 +1,78 @@
#include "FormFactors.h"
+#include "CepGen/Core/Exception.h"
+
#include "CepGen/StructureFunctions/ALLM.h"
#include "CepGen/StructureFunctions/BlockDurandHa.h"
#include "CepGen/StructureFunctions/FioreBrasse.h"
#include "CepGen/StructureFunctions/GenericLHAPDF.h"
#include "CepGen/StructureFunctions/SuriYennie.h"
#include "CepGen/StructureFunctions/SzczurekUleshchenko.h"
namespace CepGen
{
FormFactors
FormFactors::Trivial()
{
return FormFactors( 1.0, 1.0 );
}
FormFactors
FormFactors::ProtonElastic( double q2 )
{
const double mp2 = Constants::mp*Constants::mp;
const double GE = pow( 1.+q2/0.71, -2. ), GE2 = GE*GE;
const double GM = 2.79*GE, GM2 = GM*GM;
return FormFactors( ( 4.*mp2*GE2 + q2*GM2 ) / ( 4.*mp2 + q2 ), GM2 );
}
FormFactors
FormFactors::ProtonInelastic( const StructureFunctions::Type& sf, double q2, double mi2, double mf2 )
{
switch ( sf ) {
case StructureFunctions::ElasticProton:
InWarning( "Elastic proton form factors requested! Check your process definition!" );
return FormFactors::ProtonElastic( q2 );
case StructureFunctions::SuriYennie:
return FormFactors::SuriYennie( q2, mi2, mf2 );
case StructureFunctions::SzczurekUleshchenko:
return FormFactors::SzczurekUleshchenko( q2, mi2, mf2 );
case StructureFunctions::FioreBrasse:
return FormFactors::FioreBrasse( q2, mi2, mf2 );
default: throw Exception( __PRETTY_FUNCTION__, "Invalid structure functions required!", FatalError );
}
}
FormFactors
FormFactors::SuriYennie( double q2, double mi2, double mf2 )
{
const double x = q2 / ( q2 + mf2 - mi2 );
const SF::SuriYennie suriyennie, sy = suriyennie( q2, x );
//std::cout << "---> " << sy.FM << "\t" << sy.F2*x/q2 << "\t" << sy.F2*x*sqrt(mi2)/q2 << std::endl;
return FormFactors( sy.F2 * x * sqrt( mi2 ) / q2, sy.FM ); //FIXME
}
FormFactors
FormFactors::FioreBrasse( double q2, double mi2, double mf2 )
{
const double x = q2 / ( q2 + mf2 - mi2 );
SF::FioreBrasse fb, sf = fb( q2, x );
return FormFactors( sf.F2 * x / q2, -2.*sf.W1 / q2 );
}
FormFactors
FormFactors::SzczurekUleshchenko( double q2, double mi2, double mf2 )
{
const double x = q2 / ( q2 + mf2 - mi2 );
SF::SzczurekUleshchenko su, sf = su( q2, x );
return FormFactors( sf.F2 * x / q2, -2.*sf.F1 / q2 );
}
std::ostream&
operator<<( std::ostream& os, const FormFactors& ff )
{
os << Form( "Form factors: electric: Fe = %.3e ; magnetic: Fm = %.3e", ff.FE, ff.FM ).c_str();
return os;
}
}
diff --git a/CepGen/Physics/Kinematics.h b/CepGen/Physics/Kinematics.h
index af3c639..6c72ef6 100644
--- a/CepGen/Physics/Kinematics.h
+++ b/CepGen/Physics/Kinematics.h
@@ -1,100 +1,101 @@
#ifndef CepGen_Physics_Kinematics_h
#define CepGen_Physics_Kinematics_h
#include <iomanip>
#include <algorithm>
#include <string>
#include "CepGen/Core/utils.h"
+#include "CepGen/Core/Exception.h"
#include "CepGen/StructureFunctions/StructureFunctions.h"
#include "CepGen/Event/Particle.h"
#include "Cuts.h"
using std::cout;
namespace CepGen
{
/// List of kinematic constraints to apply on the process phase space.
class Kinematics
{
public:
/// Validity interval for a variable
class Limits : private std::pair<double,double>
{
public:
/// Define lower and upper limits on a quantity
Limits( double min=invalid_, double max=invalid_ ) : std::pair<double,double>( min, max ) {}
/// Lower limit to apply on the variable
double min() const { return first; }
/// Lower limit to apply on the variable
double& min() { return first; }
/// Upper limit to apply on the variable
double max() const { return second; }
/// Upper limit to apply on the variable
double& max() { return second; }
double x( double v ) const {
if ( v < 0. || v > 1. ) { InError( Form( "x must be comprised between 0 and 1 ; x value = %.5e", v ) ); }
if ( !hasMin() || !hasMax() ) return invalid_;
return first + ( second-first ) * v;
}
/// Specify the lower and upper limits on the variable
void in( double low, double up ) { first = low; second = up; }
/// Full variable range allowed
double range() const { return ( !hasMin() || !hasMax() ) ? 0. : second-first; }
/// Have a lower limit?
bool hasMin() const { return first != invalid_; }
/// Have an upper limit?
bool hasMax() const { return second != invalid_; }
/// Human-readable expression of the limits
friend std::ostream& operator<<( std::ostream&, const Limits& );
private:
static constexpr double invalid_ = -999.999;
};
public:
Kinematics();
~Kinematics();
/// Type of kinematics to consider for the process
enum ProcessMode {
ElectronProton = 0, ///< electron-proton elastic case
ElasticElastic = 1, ///< proton-proton elastic case
ElasticInelastic = 2, ///< proton-proton single-dissociative (or inelastic-elastic) case
InelasticElastic = 3, ///< proton-proton single-dissociative (or elastic-inelastic) case
InelasticInelastic = 4, ///< proton-proton double-dissociative case
ProtonElectron,
ElectronElectron
};
/// Human-readable format of a process mode (elastic/dissociative parts)
friend std::ostream& operator<<( std::ostream&, const ProcessMode& );
/// Dump all the parameters used in this process cross-section computation
/// or events generation
- void dump( std::ostream& os=std::cout ) const;
+ void dump( std::ostream& os = Logger::get().outputStream ) const;
/// Incoming particles' momentum (in \f$\text{GeV}/c\f$)
std::pair<double,double> inp;
/// Set the incoming particles' momenta (if the collision is symmetric)
inline void setSqrtS( double sqrts ) { inp = { sqrts*0.5, sqrts*0.5 }; }
/// Beam/primary particle's PDG identifier
std::pair<Particle::ParticleCode,Particle::ParticleCode> inpdg;
/// PDG id of the outgoing central particles
std::vector<Particle::ParticleCode> central_system;
/// Type of kinematics to consider for the phase space
ProcessMode mode;
/// Type of structure functions to consider
StructureFunctions::Type structure_functions;
/// Cuts on the central system produced
std::map<Cuts::Central, Limits> central_cuts;
/// Cuts on the beam remnants system
std::map<Cuts::Remnants, Limits> remnant_cuts;
/// Cuts on the initial particles kinematics
std::map<Cuts::InitialState, Limits> initial_cuts;
};
}
#endif
diff --git a/CepGen/Processes/PPtoWW.cpp b/CepGen/Processes/PPtoWW.cpp
index 120756e..4400c17 100644
--- a/CepGen/Processes/PPtoWW.cpp
+++ b/CepGen/Processes/PPtoWW.cpp
@@ -1,356 +1,356 @@
#include "PPtoWW.h"
#include <assert.h>
using namespace CepGen::Process;
PPtoWW::PPtoWW() :
- GenericKTProcess( "pptoww", "gamma,gamma->W+,W-", 4, { Particle::Photon, Particle::Photon }, { Particle::W, Particle::W } )
+ GenericKTProcess( "pptoww", "gamma,gamma->W+,W-", 4, { { Particle::Photon, Particle::Photon } }, { Particle::W, Particle::W } )
{}
void
PPtoWW::prepareKTKinematics()
{
const Kinematics::Limits rap_limits = cuts_.central_cuts[Cuts::rapidity_single];
// outgoing Ws
y1_ = rap_limits.x( xkt( 0 ) );
y2_ = rap_limits.x( xkt( 1 ) );
DebuggingInsideLoop( Form( "W bosons rapidities (%.2f < y < %.2f): %f / %f", rap_limits.min(), rap_limits.max(), y1_, y2_ ) );
Kinematics::Limits ptdiff_limits = cuts_.central_cuts[Cuts::pt_diff];
if ( !ptdiff_limits.hasMax() ) ptdiff_limits.max() = 500.; //FIXME
pt_diff_ = ptdiff_limits.x( xkt( 2 ) );
phi_pt_diff_ = 2.*M_PI*xkt( 3 );
DebuggingInsideLoop( Form( "W bosons pt difference:\n\t"
" mag = %f (%.2f < Dpt < %.2f)\n\t"
" phi = %f",
pt_diff_, ptdiff_limits.min(), ptdiff_limits.max(), phi_pt_diff_ ) );
}
double
PPtoWW::computeJacobian()
{
double jac = GenericKTProcess::minimalJacobian();
jac *= cuts_.central_cuts[Cuts::rapidity_single].range(); // d(y1)
jac *= cuts_.central_cuts[Cuts::rapidity_single].range(); // d(y2)
jac *= cuts_.central_cuts[Cuts::pt_diff].range(); // d(Dpt)
jac *= 2.*M_PI; // d(phiDpt)
return jac;
}
double
PPtoWW::computeKTFactorisedMatrixElement()
{
const double mp = Particle::massFromPDGId( Particle::Proton ), mp2 = mp*mp;
const double mw = Particle::massFromPDGId( Particle::W ), mw2 = mw*mw;
//=================================================================
// How matrix element is calculated
//=================================================================
const unsigned short method = 1;
//=================================================================
// matrix element computation
//=================================================================
//const double stild = s_/2.*(1+sqrt(1.-(4*pow(mp2, 2))/s_*s_));
// Inner photons
const double q1tx = qt1_*cos( phi_qt1_ ), q1ty = qt1_*sin( phi_qt1_ ),
q2tx = qt2_*cos( phi_qt2_ ), q2ty = qt2_*sin( phi_qt2_ );
DebuggingInsideLoop( Form( "q1t(x/y) = %e / %e\n\t"
"q2t(x/y) = %e / %e", q1tx, q1ty, q2tx, q2ty ) );
// Two-photon system
const double ptsumx = q1tx+q2tx,
ptsumy = q1ty+q2ty,
ptsum = sqrt( ptsumx*ptsumx+ptsumy*ptsumy );
const double ptdiffx = pt_diff_*cos( phi_pt_diff_ ),
ptdiffy = pt_diff_*sin( phi_pt_diff_ );
// Outgoing leptons
const double pt1x = 0.5 * ( ptsumx+ptdiffx ), pt1y = 0.5 * ( ptsumy+ptdiffy ), pt1 = sqrt( pt1x*pt1x+pt1y*pt1y ),
pt2x = 0.5 * ( ptsumx-ptdiffx ), pt2y = 0.5 * ( ptsumy-ptdiffy ), pt2 = sqrt( pt2x*pt2x+pt2y*pt2y );
const Kinematics::Limits pt_limits = cuts_.central_cuts[Cuts::pt_single];
if ( pt_limits.hasMin() && ( pt1 < pt_limits.min() || pt2 < pt_limits.min() ) ) return 0.;
if ( pt_limits.hasMax() && ( pt1 > pt_limits.max() || pt2 > pt_limits.max() ) ) return 0.;
// transverse mass for the two leptons
const double amt1 = sqrt( pt1*pt1+mw2 ),
amt2 = sqrt( pt2*pt2+mw2 );
//=================================================================
// a window in two-boson invariant mass
//=================================================================
const double invm = sqrt( amt1*amt1 + amt2*amt2 + 2.*amt1*amt2*cosh( y1_-y2_ ) - ptsum*ptsum );
const Kinematics::Limits invm_limits = cuts_.central_cuts[Cuts::mass_sum];
if ( invm_limits.hasMin() && invm < invm_limits.min() ) return 0.;
if ( invm_limits.hasMax() && invm > invm_limits.max() ) return 0.;
//=================================================================
// a window in transverse momentum difference
//=================================================================
const Kinematics::Limits ptdiff_limits = cuts_.central_cuts[Cuts::pt_diff];
if ( ptdiff_limits.hasMax() && fabs( pt1-pt2 ) > ptdiff_limits.max() ) return 0.;
//=================================================================
// a window in rapidity distance
//=================================================================
const double dely = fabs( y1_-y2_ );
const Kinematics::Limits dely_limits = cuts_.central_cuts[Cuts::dely];
if ( dely_limits.hasMin() && dely < dely_limits.min() ) return 0.;
if ( dely_limits.hasMax() && dely > dely_limits.max() ) return 0.;
//=================================================================
// auxiliary quantities
//=================================================================
const double alpha1 = amt1/sqs_*exp( y1_ ),
alpha2 = amt2/sqs_*exp( y2_ ),
beta1 = amt1/sqs_*exp( -y1_ ),
beta2 = amt2/sqs_*exp( -y2_ );
DebuggingInsideLoop( Form( "Sudakov parameters:\n\t"
" alpha1/2 = %f / %f\n\t"
" beta1/2 = %f / %f", alpha1, alpha2, beta1, beta2 ) );
const double q1t2 = q1tx*q1tx+q1ty*q1ty,
q2t2 = q2tx*q2tx+q2ty*q2ty;
//const double old_x2 = 0.; //FIXME figure out where this comes from
//const double delta_x1 = (MX_*MX_+q2t2)/((1.-old_x2)*s_);
//x1 = alpha1+alpha2+delta_x1;
const double x1 = alpha1+alpha2,
x2 = beta1 +beta2;
/*const double xi_x1 = log10(x1);
const double xi_x2 = log10(x2);*/
const double z1p = alpha1/x1, z1m = alpha2/x1,
z2p = beta1 /x2, z2m = beta2 /x2;
DebuggingInsideLoop( Form( "z(1/2)p = %f / %f\n\t"
"z(1/2)m = %f / %f", z1p, z2p, z1m, z2m ) );
if ( x1 > 1. || x2 > 1. ) return 0.; // sanity check
// FIXME FIXME FIXME
const double ak10 = event_->getOneByRole( Particle::IncomingBeam1 ).energy(),
ak1z = event_->getOneByRole( Particle::IncomingBeam1 ).momentum().pz(),
ak20 = event_->getOneByRole( Particle::IncomingBeam2 ).energy(),
ak2z = event_->getOneByRole( Particle::IncomingBeam2 ).momentum().pz();
DebuggingInsideLoop( Form( "incoming particles: p1: %f / %f\n\t"
" p2: %f / %f", ak1z, ak10, ak2z, ak20 ) );
//=================================================================
// additional conditions for energy-momentum conservation
//=================================================================
const double s1_eff = x1*s_-qt1_*qt1_, s2_eff = x2*s_-qt2_*qt2_;
DebuggingInsideLoop( Form( "s(1/2)_eff = %f / %f GeV^2\n\t"
"dilepton invariant mass = %f GeV", s1_eff, s2_eff, invm ) );
switch ( cuts_.mode ) {
case Kinematics::ElasticInelastic: if ( sqrt( s1_eff ) <= ( MY_+invm ) ) return 0.;
case Kinematics::InelasticElastic: if ( sqrt( s2_eff ) <= ( MX_+invm ) ) return 0.;
case Kinematics::InelasticInelastic: if ( sqrt( s1_eff ) <= ( MY_+invm ) ) return 0.;
if ( sqrt( s2_eff ) <= ( MX_+invm ) ) return 0.;
default: break;
}
//const double qcaptx = pcaptx, qcapty = pcapty;
//=================================================================
// four-momenta of the outgoing protons (or remnants)
//=================================================================
const double px_plus = ( 1.-x1 )*fabs( ak1z )*sqrt( 2. ),
px_minus = ( MX_*MX_ + q1t2 )*0.5/px_plus;
const double py_minus = ( 1.-x2 )*fabs( ak2z )*sqrt( 2. ), // warning! sign of pz??
py_plus = ( MY_*MY_ + q2t2 )*0.5/py_minus;
DebuggingInsideLoop( Form( "px_(+/-) = %f / %f\n\t"
"py_(+/-) = %f / %f", px_plus, px_minus, py_plus, py_minus ) );
PX_ = Particle::Momentum( -q1tx, -q1ty, ( px_plus-px_minus )/sqrt( 2. ), ( px_plus+px_minus )/sqrt( 2. ) );
PY_ = Particle::Momentum( -q2tx, -q2ty, ( py_plus-py_minus )/sqrt( 2. ), ( py_plus+py_minus )/sqrt( 2. ) );
DebuggingInsideLoop( Form( "First remnant: (E,p) = (%f, %f, %f, %f)\n\t"
"Second remnant: (E,p) = (%f, %f, %f, %f)",
PX_.px(), PX_.py(), PX_.pz(), PX_.energy(),
PY_.px(), PY_.py(), PY_.pz(), PY_.energy() ) );
/*assert( fabs( PX_.mass()-MX_ ) < 1.e-6 );
assert( fabs( PY_.mass()-MY_ ) < 1.e-6 );*/
//=================================================================
// four-momenta squared of the virtual photons
//=================================================================
const double ww = 0.5 * ( 1.+sqrt( 1.-4.*mp2/s_ ) );
// FIXME FIXME FIXME /////////////////////
Particle::Momentum q1( q1tx, q1ty, +0.5 * x1*ww*sqs_*( 1.-q1t2/x1/x1/ww/ww/s_ ), 0.5 * x1*ww*sqs_*( 1.+q1t2/x1/x1/ww/ww/s_ ) ),
q2( q2tx, q2ty, -0.5 * x2*ww*sqs_*( 1.-q2t2/x2/x2/ww/ww/s_ ), 0.5 * x2*ww*sqs_*( 1.+q2t2/x2/x2/ww/ww/s_ ) );
//////////////////////////////////////////
DebuggingInsideLoop( Form( "First photon*: (E,p), m2 = (%f, %f, %f, %f), %e\n\t"
"Second photon*: (E,p), m2 = (%f, %f, %f, %f), %e",
q1.px(), q1.py(), q1.pz(), q1.energy(), q1.mass2(),
q2.px(), q2.py(), q2.pz(), q2.energy(), q2.mass2() ) );
//const double q12 = q1.mass2(), q22 = q2.mass2();
//=================================================================
// four-momenta of the outgoing W^+ and W^-
//=================================================================
p_w1_ = Particle::Momentum( pt1x, pt1y, alpha1*ak1z + beta1*ak2z, alpha1*ak10 + beta1*ak20 );
p_w2_ = Particle::Momentum( pt2x, pt2y, alpha2*ak1z + beta2*ak2z, alpha2*ak10 + beta2*ak20 );
DebuggingInsideLoop( Form( "First W: (E,p), m = (%f, %f, %f, %f), %f\n\t"
"Second W: (E,p), m = (%f, %f, %f, %f), %f",
p_w1_.px(), p_w2_.py(), p_w1_.pz(), p_w1_.energy(), p_w1_.mass(),
p_w2_.px(), p_w2_.py(), p_w2_.pz(), p_w2_.energy(), p_w2_.mass() ) );
//assert( fabs( p_w1_.mass()-event_->getByRole( Particle::CentralSystem )[0].mass() ) < 1.e-6 );
//assert( fabs( p_w2_.mass()-event_->getByRole( Particle::CentralSystem )[1].mass() ) < 1.e-6 );
//=================================================================
// Mendelstam variables
//=================================================================
//const double shat = s_*x1*x2; // ishat = 1 (approximation)
const double shat = ( q1+q2 ).mass2(); // ishat = 2 (exact formula)
const double that1 = ( q1-p_w1_ ).mass2(), that2 = ( q2-p_w2_ ).mass2(),
uhat1 = ( q1-p_w2_ ).mass2(), uhat2 = ( q2-p_w1_ ).mass2();
DebuggingInsideLoop( Form( "that(1/2) = %f / %f\n\t"
"uhat(1/2) = %f / %f",
that1, that2, uhat1, uhat2 ) );
//const double mll = sqrt( shat );
const double that = 0.5*( that1+that2 ), uhat = 0.5*( uhat1+uhat2 );
//=================================================================
// matrix elements
//=================================================================
double amat2 = 0.;
if ( method == 0 ) {
//=================================================================
// matrix element for gamma gamma --> W^+ W^-
// (Denner+Dittmaier+Schuster)
// (work in collaboration with C. Royon)
//=================================================================
const double mw4 = mw2*mw2;
const double term1 = 2.*shat * ( 2.*shat+3.*mw2 ) / ( 3.*( mw2-that )*( mw2-uhat ) );
const double term2 = 2.*shat*shat * ( shat*shat + 3.*mw4 ) / ( 3.*( mw2-that )*( mw2-that )*( mw2-uhat )*( mw2-uhat ) );
const double auxil_gamgam = 1. - term1 + term2;
const double beta = sqrt( 1.-4.*mw2/shat );
amat2 = 3.*Constants::alphaEM*Constants::alphaEM*beta / ( 2.*shat ) * auxil_gamgam / ( beta/( 64.*M_PI*M_PI*shat ) );
}
else if ( method == 1 ) {
//=================================================================
// off-shell Nachtmann formulae
//=================================================================
const double phi_diff = phi_qt1_-phi_qt2_, phi_sum = phi_qt1_+phi_qt2_;
double amat2_0 = 0., amat2_1 = 0., amat2_interf = 0.;
for ( const auto lam3 : { -1, 0, 1 } ) {
for ( const auto lam4 : { -1, 0, 1 } ) {
double ampli_pp = WWamplitude( shat, that, +1, +1, lam3, lam4 );
double ampli_mm = WWamplitude( shat, that, -1, -1, lam3, lam4 );
double ampli_pm = WWamplitude( shat, that, +1, -1, lam3, lam4 );
double ampli_mp = WWamplitude( shat, that, -1, +1, lam3, lam4 );
amat2_0 += ampli_pp*ampli_pp + ampli_mm*ampli_mm + 2.*cos( 2.*phi_diff )*ampli_pp*ampli_mm;
amat2_1 += ampli_pm*ampli_pm + ampli_mp*ampli_mp + 2.*cos( 2.*phi_sum )*ampli_pm*ampli_mp;
amat2_interf += -2.*( cos( phi_sum+phi_diff )*( ampli_pp*ampli_pm+ampli_mm*ampli_mp ) + cos( phi_sum-phi_diff )*( ampli_pp*ampli_mp+ampli_mm*ampli_pm ) );
}
}
amat2 = amat2_0 + amat2_1 + amat2_interf;
}
//============================================
// unintegrated photon distributions
//============================================
GenericKTProcess::computeIncomingFluxes( x1, q1t2, x2, q2t2 );
//=================================================================
// factor 2.*pi from integration over phi_sum
// factor 1/4 from jacobian of transformations
// factors 1/pi and 1/pi due to integration over
// d^2 kappa_1 d^2 kappa_2 instead d kappa_1^2 d kappa_2^2
//=================================================================
const double aintegral = amat2 * ( 2.*M_PI )/ ( 16.*M_PI*M_PI*( x1*x2*s_ )*( x1*x2*s_ ) )
* flux1_/M_PI * flux2_/M_PI * 0.25
* Constants::GeV2toBarn * 0.5 / M_PI;
/*const double aintegral = amat2 / ( 16.*M_PI*M_PI*x1*x1*x2*x2*s_*s_ )
* flux1_/M_PI * flux2_/M_PI
* Constants::GeV2toBarn * 0.25;*/
//=================================================================
return aintegral*qt1_*qt2_*pt_diff_;
//=================================================================
}
void
PPtoWW::fillCentralParticlesKinematics()
{
// randomise the charge of the outgoing leptons
int sign = ( drand()>.5 ) ? +1 : -1;
//=================================================================
// first outgoing lepton
//=================================================================
Particle& ow1 = event_->getByRole( Particle::CentralSystem )[0];
ow1.setPdgId( ow1.pdgId(), sign );
ow1.setStatus( Particle::FinalState );
ow1.setMomentum( p_w1_ );
//=================================================================
// second outgoing lepton
//=================================================================
Particle& ow2 = event_->getByRole( Particle::CentralSystem )[1];
ow2.setPdgId( ow2.pdgId(), -sign );
ow2.setStatus( Particle::FinalState );
ow2.setMomentum( p_w2_ );
}
double
PPtoWW::WWamplitude( double shat, double that, short lam1, short lam2, short lam3, short lam4 ) const
{
const double mw = Particle::massFromPDGId( Particle::W ), mw2 = mw*mw;
const double cos_theta = ( 1.+2.*( that-mw2 )/shat ) / sqrt( 1. + 1.e-10 - 4.*mw2/shat );
const double sin_theta = sqrt( 1.-cos_theta*cos_theta );
const double beta = sqrt( 1.-4.*mw2/shat );
const double gam = 1./sqrt( 1.-beta*beta );
const double A = ( 1.-beta*beta*cos_theta*cos_theta );
const double e2 = 4.*M_PI*Constants::alphaEM;
const double term1 = 1./( gam*gam*A )*( ( gam*gam+1. )*( 1.-lam1*lam2 )* sin_theta*sin_theta - ( 1.+lam1*lam2 ) );
const double term2 = -sqrt( 2. )/( gam*A )*( lam1-lam2 ) * ( 1.+lam1*lam3*cos_theta )*sin_theta;
const double term3 = -0.5/A*( 2.*beta*( lam1+lam2 )*( lam3+lam4 ) - ( 1./gam/gam )*( 1.+lam3*lam4 )*( 2.*lam1*lam2+( 1.-lam1*lam2 ) * cos_theta*cos_theta ) +( 1.+lam1*lam2*lam3*lam4 )*( 3.+lam1*lam2 ) + 2.*( lam1-lam2 )*( lam3-lam4 )*cos_theta + ( 1.-lam1*lam2 )*( 1.-lam3*lam4 )*cos_theta*cos_theta );
const double term4 = -sqrt( 2. )/( gam*A )*( lam2-lam1 ) *( 1.+lam2*lam4*cos_theta )*sin_theta;
if ( lam3 == 0 && lam4 == 0 ) return e2*term1;
if ( lam4 == 0 ) return e2*term2;
if ( lam3 == 0 ) return e2*term4;
if ( lam3 != 0 && lam4 != 0 ) return e2*term3;
return 0.;
}
diff --git a/CepGen/StructureFunctions/ChristyBosted.cpp b/CepGen/StructureFunctions/ChristyBosted.cpp
index 54134d7..f4e0eb1 100644
--- a/CepGen/StructureFunctions/ChristyBosted.cpp
+++ b/CepGen/StructureFunctions/ChristyBosted.cpp
@@ -1,259 +1,260 @@
#include "ChristyBosted.h"
#include "CepGen/Event/Particle.h"
+#include "CepGen/Core/Exception.h"
namespace CepGen
{
namespace SF
{
double
ChristyBosted::resmod507( char sf, double w2, double q2 ) const
{
const double mp = Constants::mp, mp2 = mp*mp,
mpi = Constants::mpi, mpi2 = mpi*mpi,
meta = Particle::massFromPDGId( Particle::Eta ), meta2 = meta*meta;
const double w = sqrt( w2 );
const double xb = q2/( q2+w2-mp2 );
double m0 = 0., q20 = 0.;
if ( sf == 'T' ) { // transverse
m0 = 0.125;
q20 = 0.05;
}
else if ( sf == 'L' ) {
m0 = params_.m0;
q20 = 0.125;
}
else {
InError( "Invalid direction retrieved! Aborting." )
return 0.;
}
const double norm_q2 = 1./0.330/0.330;
const double t = log( log( ( q2+m0 )*norm_q2 )/log( m0*norm_q2 ) );
//--- calculate kinematics needed for threshold relativistic B-W
// equivalent photon energies
const double k = 0.5 * ( w2 - mp2 )/mp;
const double kcm = 0.5 * ( w2 - mp2 )/w;
const double epicm = 0.5 * ( w2 + mpi2 - mp2 )/w, ppicm = sqrt( std::max( 0., epicm* epicm - mpi2 ) );
const double epi2cm = 0.5 * ( w2 + 4.*mpi2 - mp2 )/w, ppi2cm = sqrt( std::max( 0., epi2cm*epi2cm - 4*mpi2 ) );
const double eetacm = 0.5 * ( w2 + meta2 - mp2 )/w, petacm = sqrt( std::max( 0., eetacm*eetacm - meta2 ) );
std::array<double,7> width, height, pgam;
for ( unsigned short i = 0; i < 7; ++i ) {
const Parameterisation::ResonanceParameters& res = params_.resonances[i];
width[i] = res.width;
//--- calculate partial widths
//----- 1-pion decay mode
const double x02 = res.x0*res.x0;
const double partial_width_singlepi = pow( ppicm /res.pcmr( mpi2 ), 2.*res.angular_momentum+1. )
* pow( ( res.pcmr( mpi2 )*res.pcmr( mpi2 )+x02 )/( ppicm *ppicm +x02 ), res.angular_momentum );
//----- 2-pion decay mode
const double partial_width_doublepi = pow( ppi2cm/res.pcmr( 4.*mpi2 ), 2.*( res.angular_momentum+2. ) )
* pow( ( res.pcmr( 4.*mpi2 )*res.pcmr( 4.*mpi2 )+x02 )/( ppi2cm*ppi2cm+x02 ), res.angular_momentum+2 )
* w / res.mass;
//----- eta decay mode (only for S11's)
const double partial_width_eta = ( res.br.eta == 0. ) ? 0. :
pow( petacm/res.pcmr( meta2 ), 2.*res.angular_momentum+1. )
* pow( ( res.pcmr( meta2 )*res.pcmr( meta2 )+x02 )/( petacm*petacm+x02 ), res.angular_momentum );
// virtual photon width
pgam[i] = res.width * pow( kcm/res.kcmr(), 2 ) * ( res.kcmr()*res.kcmr()+x02 )/( kcm*kcm+x02 );
width[i] = ( partial_width_singlepi * res.br.singlepi
+ partial_width_doublepi * res.br.doublepi
+ partial_width_eta * res.br.eta ) * res.width;
//--- resonance Q^2 dependence calculations
if ( sf == 'T' ) height[i] = res.A0_T*( 1.+res.fit_parameters[0]*q2/( 1.+res.fit_parameters[1]*q2 ) )/pow( 1.+q2/0.91, res.fit_parameters[2] );
else if ( sf == 'L' ) height[i] = res.A0_L/( 1.+res.fit_parameters[3]*q2 )*q2*exp( -q2*res.fit_parameters[4] );
height[i] = height[i]*height[i];
}
//--- calculate Breit-Wigners for all resonances
double sig_res = 0.;
for ( unsigned short i = 0; i < 7; ++i ) {
const Parameterisation::ResonanceParameters res = params_.resonances[i];
const double mass2 = res.mass*res.mass, width2 = width[i]*width[i];
const double sigr = height[i]*res.kr()/k*res.kcmr()/kcm/res.width * ( width[i]*pgam[i] / ( pow( w2-mass2, 2 ) + mass2*width2 ) );
sig_res += sigr;
}
sig_res *= w;
//--- non-resonant background calculation
const double xpr = 1./( 1.+( w2-pow( mp+mpi, 2 ) )/( q2+q20 ) );
if ( xpr > 1. ) return 0.; // FIXME
double sig_nr = 0.;
if ( sf == 'T' ) { // transverse
const double wdif = w - ( mp + mpi );
if ( wdif >= 0. ) {
for ( unsigned short i = 0; i < 2; ++i ) {
const double expo = params_.continuum.transverse[i].fit_parameters[1]
+ params_.continuum.transverse[i].fit_parameters[2]*q2
+ params_.continuum.transverse[i].fit_parameters[3]*q2*q2;
sig_nr += params_.continuum.transverse[i].sig0 / pow( q2+params_.continuum.transverse[i].fit_parameters[0], expo ) * pow( wdif, i+1.5 );
}
}
sig_nr *= xpr;
}
else if ( sf == 'L' ) { // longitudinal
for ( unsigned short i = 0; i < 1; ++i ) {
const double expo = params_.continuum.longitudinal[i].fit_parameters[0]
+ params_.continuum.longitudinal[i].fit_parameters[1];
sig_nr += params_.continuum.longitudinal[i].sig0
* pow( 1.-xpr, expo )/( 1.-xb )
* pow( q2/( q2+q20 ), params_.continuum.longitudinal[i].fit_parameters[2] )/( q2+q20 )
* pow( xpr, params_.continuum.longitudinal[i].fit_parameters[3]+params_.continuum.longitudinal[i].fit_parameters[4]*t );
}
}
return sig_res + sig_nr;
}
ChristyBosted::Parameterisation
ChristyBosted::Parameterisation::standard()
{
Parameterisation params;
params.m0 = 4.2802;
params.continuum.transverse = { {
ContinuumParameters::DirectionParameters( 246.06, { { 0.067469, 1.3501, 0.12054, -0.0038495 } } ),
ContinuumParameters::DirectionParameters( -89.360, { { 0.20977, 1.5715, 0.090736, 0.010362 } } )
} };
params.continuum.longitudinal = { {
ContinuumParameters::DirectionParameters( 86.746, { { 0., 4.0294, 3.1285, 0.33403, 4.9623 } } )
} };
//--- P33(1232)
ResonanceParameters p33;
p33.br = ResonanceParameters::BranchingRatios( 1., 0., 0. );
p33.angular_momentum = 1.;
//p33.x0 = 0.15;
p33.x0 = 0.14462;
p33.mass = 1.2298;
p33.width = 0.13573;
p33.fit_parameters = { { 4.2291, 1.2598, 2.1242, 19.910, 0.22587 } };
p33.A0_T = 7.7805;
p33.A0_L = 29.414;
params.resonances.emplace_back( p33 );
//--- S11(1535)
ResonanceParameters s11_1535;
s11_1535.br = ResonanceParameters::BranchingRatios( 0.45, 0.1, 0.45 );
s11_1535.angular_momentum = 0.;
s11_1535.x0 = 0.215;
s11_1535.mass = 1.5304;
s11_1535.width = 0.220;
s11_1535.fit_parameters = { { 6823.2, 33521., 2.5686, 0., 0. } };
s11_1535.A0_T = 6.3351;
s11_1535.A0_L = 0.;
params.resonances.emplace_back( s11_1535 );
//--- D13(1520)
ResonanceParameters d13;
d13.br = ResonanceParameters::BranchingRatios( 0.65, 0.35, 0. );
d13.angular_momentum = 2.;
d13.x0 = 0.215;
d13.mass = 1.5057;
d13.width = 0.082956;
d13.fit_parameters = { { 21.240, 0.055746, 2.4886, 97.046, 0.31042 } };
d13.A0_T = 0.60347;
d13.A0_L = 157.92;
params.resonances.emplace_back( d13 );
//--- F15(1680)
ResonanceParameters f15;
f15.br = ResonanceParameters::BranchingRatios( 0.65, 0.35, 0. );
f15.angular_momentum = 3.;
f15.x0 = 0.215;
f15.mass = 1.6980;
f15.width = 0.095782;
f15.fit_parameters = { { -0.28789, 0.18607, 0.063534, 0.038200, 1.2182 } };
f15.A0_T = 2.3305;
f15.A0_L = 4.2160;
params.resonances.emplace_back( f15 );
//--- S11(1650)
ResonanceParameters s11_1650;
s11_1650.br = ResonanceParameters::BranchingRatios( 0.4, 0.5, 0.1 );
s11_1650.angular_momentum = 0.;
s11_1650.x0 = 0.215;
s11_1650.mass = 1.6650;
s11_1650.width = 0.10936;
s11_1650.fit_parameters = { { -0.56175, 0.38964, 0.54883, 0.31393, 2.9997 } };
s11_1650.A0_T = 1.9790;
s11_1650.A0_L = 13.764;
params.resonances.emplace_back( s11_1650 );
//--- P11(1440) roper
ResonanceParameters p11;
p11.br = ResonanceParameters::BranchingRatios( 0.65, 0.35, 0. );
p11.angular_momentum = 1.;
p11.x0 = 0.215;
p11.mass = 1.4333;
p11.width = 0.37944;
p11.fit_parameters = { { 46.213, 0.19221, 1.9141, 0.053743, 1.3091 } };
p11.A0_T = 0.022506;
p11.A0_L = 5.5124;
params.resonances.emplace_back( p11 );
//--- F37(1950)
ResonanceParameters f37;
f37.br = ResonanceParameters::BranchingRatios( 0.5, 0.5, 0. );
f37.angular_momentum = 3.;
f37.x0 = 0.215;
f37.mass = 1.9341;
f37.width = 0.380;
f37.fit_parameters = { { 0., 0., 1., 1.8951, 0.51376 } };
f37.A0_T = 3.4187;
f37.A0_L = 1.8951;
params.resonances.emplace_back( f37 );
return params;
}
ChristyBosted
ChristyBosted::operator()( double q2, double xbj ) const
{
const double mp2 = Constants::mp*Constants::mp;
const double w2 = mp2 + q2*( 1.-xbj )/xbj;
const double w_min = Constants::mp+Constants::mpi;
ChristyBosted cb;
if ( sqrt( w2 ) < w_min ) return cb;
//-----------------------------
// modification of Christy-Bosted at large q2 as described in the LUXqed paper
//-----------------------------
const double q21 = 30., q20 = 8.;
const double delq2 = q2 - q20;
const double qq = q21 - q20;
const double prefac = 1./( 4.*M_PI*M_PI*Constants::alphaEM ) * ( 1.-xbj );
//------------------------------
double q2_eff = q2, w2_eff = w2;
if ( q2 > q20 ) {
q2_eff = q20 + delq2/( 1.+delq2/qq );
w2_eff = mp2 + q2_eff*( 1.-xbj )/xbj;
}
const double tau = 4.*xbj*xbj*mp2/q2_eff;
const double sigT = resmod507( 'T', w2_eff, q2_eff ), sigL = resmod507( 'L', w2_eff, q2_eff );
cb.F2 = prefac * q2_eff / ( 1+tau ) * ( sigT+sigL ) / Constants::GeV2toBarn * 1.e6;
if ( q2 > q20 ) cb.F2 *= q21/( q21 + delq2 );
if ( sigT != 0. ) {
const double R = sigL/sigT;
cb.FL = cb.F2*( 1+tau )*R/( 1.+R );
}
return cb;
}
}
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 8:23 PM (1 d, 5 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3806022
Default Alt Text
(125 KB)

Event Timeline