diff --git a/CepGen/Cards/LpairHandler.cpp b/CepGen/Cards/LpairHandler.cpp
index 1682054..fc81d0c 100644
--- a/CepGen/Cards/LpairHandler.cpp
+++ b/CepGen/Cards/LpairHandler.cpp
@@ -1,217 +1,217 @@
#include "CepGen/Cards/LpairHandler.h"
#include "CepGen/Core/ParametersList.h"
#include "CepGen/Core/Exception.h"
#include "CepGen/Physics/PDG.h"
#include "CepGen/StructureFunctions/StructureFunctions.h"
#include "CepGen/StructureFunctions/LHAPDF.h"
#include "CepGen/Processes/ProcessesHandler.h"
#include "CepGen/Hadronisers/Pythia8Hadroniser.h"
#include <fstream>
namespace CepGen
{
namespace Cards
{
const int LpairHandler::kInvalid = 99999;
//----- specialization for LPAIR input cards
LpairHandler::LpairHandler( const char* file ) :
proc_params_( new ParametersList ),
str_fun_( 11 ), hi_1_( { 0, 0 } ), hi_2_( { 0, 0 } )
{
std::ifstream f( file, std::fstream::in );
if ( !f.is_open() )
throw CG_FATAL( "LpairHandler" ) << "Failed to parse file \"" << file << "%s\".";
init( ¶ms_ );
//--- parse all fields
std::unordered_map<std::string, std::string> m_params;
std::string key, value;
std::ostringstream os;
while ( f >> key >> value ) {
if ( key[0] == '#' ) // FIXME need to ensure there is no extra space before!
continue;
setParameter( key, value );
m_params.insert( { key, value } );
if ( getDescription( key ) != "null" )
os << "\n>> " << key << " = " << std::setw( 15 ) << getParameter( key )
<< " (" << getDescription( key ) << ")";
}
f.close();
//--- parse the process name
auto proc = CepGen::ProcessesHandler::get().build( proc_name_, *proc_params_ );
params_.setProcess( std::move( proc ) );
//--- parse the structure functions code
const unsigned long kLHAPDFCodeDec = 10000000, kLHAPDFPartDec = 1000000;
if ( str_fun_ / kLHAPDFCodeDec == 1 ) { // SF from parton
params_.kinematics.structure_functions = sf::Parameterisation::build( sf::Type::LHAPDF );
auto sf = dynamic_cast<sf::LHAPDF*>( params_.kinematics.structure_functions.get() );
const unsigned long icode = str_fun_ % kLHAPDFCodeDec;
sf->params.pdf_code = icode % kLHAPDFPartDec;
sf->params.mode = (sf::LHAPDF::Parameters::Mode)( icode / kLHAPDFPartDec ); // 0, 1, 2
}
else
params_.kinematics.structure_functions = sf::Parameterisation::build( (sf::Type)str_fun_ );
//--- parse the integration algorithm name
if ( integr_type_ == "plain" )
params_.integrator.type = Integrator::Type::plain;
else if ( integr_type_ == "Vegas" )
params_.integrator.type = Integrator::Type::Vegas;
else if ( integr_type_ == "MISER" )
params_.integrator.type = Integrator::Type::MISER;
else if ( integr_type_ != "" )
throw CG_FATAL( "LpairHandler" ) << "Unrecognized integrator type: " << integr_type_ << "!";
//--- parse the hadronisation algorithm name
if ( hadr_name_ == "pythia8" )
- params_.setHadroniser( new Hadroniser::Pythia8Hadroniser( params_, ParametersList() ) );
+ params_.setHadroniser( new hadroniser::Pythia8Hadroniser( params_, ParametersList() ) );
if ( m_params.count( "IEND" ) )
setValue<bool>( "IEND", ( std::stoi( m_params["IEND"] ) > 1 ) );
//--- check if we are dealing with heavy ions for incoming states
HeavyIon hi1{ hi_1_.first, (Element)hi_1_.second }, hi2{ hi_2_.first, (Element)hi_2_.second };
if ( hi1 )
params_.kinematics.incoming_beams.first.pdg = hi1;
if ( hi2 )
params_.kinematics.incoming_beams.second.pdg = hi2;
CG_INFO( "LpairHandler" ) << "File '" << file << "' succesfully opened!\n\t"
<< "The following parameters are set:" << os.str();
}
void
LpairHandler::init( Parameters* params )
{
//-------------------------------------------------------------------------------------------
// Process/integration/hadronisation parameters
//-------------------------------------------------------------------------------------------
registerParameter<std::string>( "PROC", "Process name to simulate", &proc_name_ );
registerParameter<std::string>( "ITYP", "Integration algorithm", &integr_type_ );
registerParameter<std::string>( "HADR", "Hadronisation algorithm", &hadr_name_ );
registerParameter<std::string>( "KMRG", "KMR grid interpolation path", ¶ms_.kinematics.kmr_grid_path );
//-------------------------------------------------------------------------------------------
// General parameters
//-------------------------------------------------------------------------------------------
registerParameter<bool>( "IEND", "Generation type", ¶ms->generation.enabled );
registerParameter<bool>( "NTRT", "Smoothen the integrand", ¶ms->generation.treat );
registerParameter<int>( "DEBG", "Debugging verbosity", (int*)&Logger::get().level );
registerParameter<int>( "NCVG", "Number of function calls", (int*)¶ms->integrator.ncvg );
registerParameter<int>( "ITVG", "Number of integration iterations", (int*)¶ms->integrator.vegas.iterations );
registerParameter<int>( "SEED", "Random generator seed", (int*)¶ms->integrator.rng_seed );
registerParameter<int>( "NTHR", "Number of threads to use for events generation", (int*)¶ms->generation.num_threads );
registerParameter<int>( "MODE", "Subprocess' mode", (int*)¶ms->kinematics.mode );
registerParameter<int>( "NCSG", "Number of points to probe", (int*)¶ms->generation.num_points );
registerParameter<int>( "NGEN", "Number of events to generate", (int*)¶ms->generation.maxgen );
registerParameter<int>( "NPRN", "Number of events before printout", (int*)¶ms->generation.gen_print_every );
//-------------------------------------------------------------------------------------------
// Process-specific parameters
//-------------------------------------------------------------------------------------------
registerParameter<int>( "METH", "Computation method (kT-factorisation)", &proc_params_->operator[]<int>( "method" ) );
registerParameter<int>( "IPOL", "Polarisation states to consider", &proc_params_->operator[]<int>( "polarisationStates" ) );
//-------------------------------------------------------------------------------------------
// Process kinematics parameters
//-------------------------------------------------------------------------------------------
registerParameter<int>( "PMOD", "Outgoing primary particles' mode", &str_fun_ );
registerParameter<int>( "EMOD", "Outgoing primary particles' mode", &str_fun_ );
registerParameter<int>( "PAIR", "Outgoing particles' PDG id", (int*)&proc_params_->operator[]<int>( "pair" ) );
registerParameter<int>( "INA1", "Heavy ion atomic weight (1st incoming beam)", (int*)&hi_1_.first );
registerParameter<int>( "INZ1", "Heavy ion atomic number (1st incoming beam)", (int*)&hi_1_.second );
registerParameter<int>( "INA2", "Heavy ion atomic weight (1st incoming beam)", (int*)&hi_2_.first );
registerParameter<int>( "INZ2", "Heavy ion atomic number (1st incoming beam)", (int*)&hi_2_.second );
registerParameter<double>( "INP1", "Momentum (1st primary particle)", ¶ms->kinematics.incoming_beams.first.pz );
registerParameter<double>( "INP2", "Momentum (2nd primary particle)", ¶ms->kinematics.incoming_beams.second.pz );
registerParameter<double>( "INPP", "Momentum (1st primary particle)", ¶ms->kinematics.incoming_beams.first.pz );
registerParameter<double>( "INPE", "Momentum (2nd primary particle)", ¶ms->kinematics.incoming_beams.second.pz );
registerParameter<double>( "PTCT", "Minimal transverse momentum (single central outgoing particle)", ¶ms->kinematics.cuts.central.pt_single.min() );
registerParameter<double>( "MSCT", "Minimal central system mass", ¶ms->kinematics.cuts.central.mass_sum.min() );
registerParameter<double>( "ECUT", "Minimal energy (single central outgoing particle)", ¶ms->kinematics.cuts.central.energy_single.min() );
registerParameter<double>( "ETMN", "Minimal pseudo-rapidity (central outgoing particles)", ¶ms->kinematics.cuts.central.eta_single.min() );
registerParameter<double>( "ETMX", "Maximal pseudo-rapidity (central outgoing particles)", ¶ms->kinematics.cuts.central.eta_single.max() );
registerParameter<double>( "YMIN", "Minimal rapidity (central outgoing particles)", ¶ms->kinematics.cuts.central.rapidity_single.min() );
registerParameter<double>( "YMAX", "Maximal rapidity (central outgoing particles)", ¶ms->kinematics.cuts.central.rapidity_single.max() );
registerParameter<double>( "Q2MN", "Minimal Q² = -q² (exchanged parton)", ¶ms->kinematics.cuts.initial.q2.min() );
registerParameter<double>( "Q2MX", "Maximal Q² = -q² (exchanged parton)", ¶ms->kinematics.cuts.initial.q2.max() );
registerParameter<double>( "MXMN", "Minimal invariant mass of proton remnants", ¶ms->kinematics.cuts.remnants.mass_single.min() );
registerParameter<double>( "MXMX", "Maximal invariant mass of proton remnants", ¶ms->kinematics.cuts.remnants.mass_single.max() );
}
void
LpairHandler::store( const char* file )
{
std::ofstream f( file, std::fstream::out | std::fstream::trunc );
if ( !f.is_open() ) {
CG_ERROR( "LpairHandler" ) << "Failed to open file \"" << file << "%s\" for writing.";
}
for ( const auto& it : p_strings_ )
if ( it.second.value )
f << it.first << " = " << *it.second.value << "\n";
for ( const auto& it : p_ints_ )
if ( it.second.value )
f << it.first << " = " << *it.second.value << "\n";
for ( const auto& it : p_doubles_ )
if ( it.second.value )
f << it.first << " = " << *it.second.value << "\n";
for ( const auto& it : p_bools_ )
if ( it.second.value )
f << it.first << " = " << *it.second.value << "\n";
f.close();
}
void
LpairHandler::setParameter( const std::string& key, const std::string& value )
{
try { setValue<double>( key.c_str(), std::stod( value ) ); } catch ( std::invalid_argument& ) {}
try { setValue<int>( key.c_str(), std::stoi( value ) ); } catch ( std::invalid_argument& ) {}
//setValue<bool>( key.c_str(), std::stoi( value ) );
setValue<std::string>( key.c_str(), value );
}
std::string
LpairHandler::getParameter( std::string key ) const
{
double dd = getValue<double>( key.c_str() );
if ( dd != -999. )
return std::to_string( dd );
int ui = getValue<int>( key.c_str() );
if ( ui != 999 )
return std::to_string( ui );
//if ( out = getValue<bool>( key.c_str() ) );
return getValue<std::string>( key.c_str() );
}
std::string
LpairHandler::getDescription( std::string key ) const
{
if ( p_strings_.count( key ) )
return p_strings_.find( key )->second.description;
if ( p_ints_.count( key ) )
return p_ints_.find( key )->second.description;
if ( p_doubles_.count( key ) )
return p_doubles_.find( key )->second.description;
if ( p_bools_.count( key ) )
return p_bools_.find( key )->second.description;
return "null";
}
}
}
diff --git a/CepGen/Cards/PythonHandler.cpp b/CepGen/Cards/PythonHandler.cpp
index 3393003..83a7f77 100644
--- a/CepGen/Cards/PythonHandler.cpp
+++ b/CepGen/Cards/PythonHandler.cpp
@@ -1,404 +1,404 @@
#include "CepGen/Cards/PythonHandler.h"
#include "CepGen/Core/Exception.h"
#ifdef PYTHON
#include "CepGen/Core/TamingFunction.h"
#include "CepGen/Core/Exception.h"
#include "CepGen/Core/ParametersList.h"
#include "CepGen/Processes/ProcessesHandler.h"
#include "CepGen/StructureFunctions/StructureFunctions.h"
#include "CepGen/StructureFunctions/LHAPDF.h"
#include "CepGen/StructureFunctions/MSTWGrid.h"
#include "CepGen/StructureFunctions/Schaefer.h"
#include "CepGen/Hadronisers/Pythia8Hadroniser.h"
#include <algorithm>
#if PY_MAJOR_VERSION < 3
# define PYTHON2
#endif
namespace CepGen
{
namespace Cards
{
//----- specialization for CepGen input cards
PythonHandler::PythonHandler( const char* file )
{
setenv( "PYTHONPATH", ".:..:Cards", 1 );
std::string filename = getPythonPath( file );
const size_t fn_len = filename.length()+1;
//Py_DebugFlag = 1;
//Py_VerboseFlag = 1;
#ifdef PYTHON2
char* sfilename = new char[fn_len];
snprintf( sfilename, fn_len, "%s", filename.c_str() );
#else
wchar_t* sfilename = new wchar_t[fn_len];
swprintf( sfilename, fn_len, L"%s", filename.c_str() );
#endif
if ( sfilename )
Py_SetProgramName( sfilename );
Py_InitializeEx( 1 );
if ( sfilename )
delete [] sfilename;
if ( !Py_IsInitialized() )
throw CG_FATAL( "PythonHandler" ) << "Failed to initialise the Python cards parser!";
CG_INFO( "PythonHandler" )
<< "Initialised the Python cards parser\n\t"
<< "Python version: " << Py_GetVersion() << "\n\t"
<< "Platform: " << Py_GetPlatform() << ".";
PyObject* cfg = PyImport_ImportModule( filename.c_str() ); // new
if ( !cfg )
throwPythonError( Form( "Failed to parse the configuration card %s", file ) );
PyObject* process = PyObject_GetAttrString( cfg, PROCESS_NAME ); // new
if ( !process )
throwPythonError( Form( "Failed to extract a \"%s\" keyword from the configuration card %s", PROCESS_NAME, file ) );
//--- list of process-specific parameters
ParametersList proc_params;
fillParameter( process, "processParameters", proc_params );
//--- type of process to consider
PyObject* pproc_name = getElement( process, MODULE_NAME ); // borrowed
if ( !pproc_name )
throwPythonError( Form( "Failed to extract the process name from the configuration card %s", file ) );
const std::string proc_name = get<std::string>( pproc_name );
//--- process mode
params_.kinematics.mode = (KinematicsMode)proc_params.get<int>( "mode", (int)KinematicsMode::invalid );
auto proc = CepGen::ProcessesHandler::get().build( proc_name, proc_params );
params_.setProcess( std::move( proc ) );
//--- process kinematics
PyObject* pin_kinematics = getElement( process, "inKinematics" ); // borrowed
if ( pin_kinematics )
parseIncomingKinematics( pin_kinematics );
PyObject* pout_kinematics = getElement( process, "outKinematics" ); // borrowed
if ( pout_kinematics )
parseOutgoingKinematics( pout_kinematics );
//--- taming functions
PyObject* ptam = getElement( process, "tamingFunctions" ); // borrowed
if ( ptam )
parseTamingFunctions( ptam );
Py_CLEAR( process );
PyObject* plog = PyObject_GetAttrString( cfg, "logger" ); // new
if ( plog ) {
parseLogging( plog );
Py_CLEAR( plog );
}
//--- hadroniser parameters
PyObject* phad = PyObject_GetAttrString( cfg, "hadroniser" ); // new
if ( phad ) {
parseHadroniser( phad );
Py_CLEAR( phad );
}
//--- generation parameters
PyObject* pint = PyObject_GetAttrString( cfg, "integrator" ); // new
if ( pint ) {
parseIntegrator( pint );
Py_CLEAR( pint );
}
PyObject* pgen = PyObject_GetAttrString( cfg, "generator" ); // new
if ( pgen ) {
parseGenerator( pgen );
Py_CLEAR( pgen );
}
//--- finalisation
Py_CLEAR( cfg );
}
PythonHandler::~PythonHandler()
{
if ( Py_IsInitialized() )
Py_Finalize();
}
void
PythonHandler::parseIncomingKinematics( PyObject* kin )
{
//--- retrieve the beams PDG ids
std::vector<double> beams_pz;
fillParameter( kin, "pz", beams_pz );
if ( beams_pz.size() == 2 ) {
params_.kinematics.incoming_beams.first.pz = beams_pz.at( 0 );
params_.kinematics.incoming_beams.second.pz = beams_pz.at( 1 );
}
//--- retrieve the beams longitudinal momentum
std::vector<int> beams_pdg;
fillParameter( kin, "pdgIds", beams_pdg );
if ( beams_pdg.size() == 2 ) {
params_.kinematics.incoming_beams.first.pdg = (PDG)beams_pdg.at( 0 );
params_.kinematics.incoming_beams.second.pdg = (PDG)beams_pdg.at( 1 );
}
double sqrt_s = -1.;
fillParameter( kin, "cmEnergy", sqrt_s );
fillParameter( kin, "kmrGridPath", params_.kinematics.kmr_grid_path );
if ( sqrt_s != -1. )
params_.kinematics.setSqrtS( sqrt_s );
PyObject* psf = getElement( kin, "structureFunctions" ); // borrowed
if ( psf )
parseStructureFunctions( psf, params_.kinematics.structure_functions );
std::vector<int> kt_fluxes;
fillParameter( kin, "ktFluxes", kt_fluxes );
if ( kt_fluxes.size() > 0 )
params_.kinematics.incoming_beams.first.kt_flux = (KTFlux)kt_fluxes.at( 0 );
if ( kt_fluxes.size() > 1 )
params_.kinematics.incoming_beams.second.kt_flux = (KTFlux)kt_fluxes.at( 1 );
std::vector<int> hi_beam1, hi_beam2;
fillParameter( kin, "heavyIonA", hi_beam1 );
if ( hi_beam1.size() == 2 )
params_.kinematics.incoming_beams.first.pdg = HeavyIon{ (unsigned short)hi_beam1[0], (Element)hi_beam1[1] };
fillParameter( kin, "heavyIonB", hi_beam2 );
if ( hi_beam2.size() == 2 )
params_.kinematics.incoming_beams.second.pdg = HeavyIon{ (unsigned short)hi_beam2[0], (Element)hi_beam2[1] };
}
void
PythonHandler::parseStructureFunctions( PyObject* psf, std::shared_ptr<sf::Parameterisation>& sf_handler )
{
int str_fun = 0;
fillParameter( psf, "id", str_fun );
sf_handler = sf::Parameterisation::build( (sf::Type)str_fun );
switch( (sf::Type)str_fun ) {
case sf::Type::LHAPDF: {
auto sf = std::dynamic_pointer_cast<sf::LHAPDF>( params_.kinematics.structure_functions );
fillParameter( psf, "pdfSet", sf->params.pdf_set );
fillParameter( psf, "numFlavours", (unsigned int&)sf->params.num_flavours );
fillParameter( psf, "pdfMember", (unsigned int&)sf->params.pdf_member );
fillParameter( psf, "mode", (unsigned int&)sf->params.mode );
} break;
case sf::Type::MSTWgrid: {
auto sf = std::dynamic_pointer_cast<mstw::Grid>( params_.kinematics.structure_functions );
fillParameter( psf, "gridPath", sf->params.grid_path );
} break;
case sf::Type::Schaefer: {
auto sf = std::dynamic_pointer_cast<sf::Schaefer>( params_.kinematics.structure_functions );
fillParameter( psf, "Q2cut", sf->params.q2_cut );
std::vector<double> w2_lims;
fillParameter( psf, "W2limits", w2_lims );
if ( w2_lims.size() != 0 ) {
if ( w2_lims.size() != 2 )
throwPythonError( Form( "Invalid size for W2limits attribute: %d != 2!", w2_lims.size() ) );
else {
sf->params.w2_lo = *std::min_element( w2_lims.begin(), w2_lims.end() );
sf->params.w2_hi = *std::max_element( w2_lims.begin(), w2_lims.end() );
}
}
PyObject* pcsf = getElement( psf, "continuumSF" ); // borrowed
if ( pcsf )
parseStructureFunctions( pcsf, sf->params.continuum_model );
PyObject* ppsf = getElement( psf, "perturbativeSF" ); // borrowed
if ( ppsf )
parseStructureFunctions( ppsf, sf->params.perturbative_model );
PyObject* prsf = getElement( psf, "resonancesSF" ); // borrowed
if ( prsf )
parseStructureFunctions( prsf, sf->params.resonances_model );
fillParameter( psf, "higherTwist", (bool&)sf->params.higher_twist );
} break;
default: break;
}
}
void
PythonHandler::parseOutgoingKinematics( PyObject* kin )
{
PyObject* pparts = getElement( kin, "minFinalState" ); // borrowed
if ( pparts && PyTuple_Check( pparts ) )
for ( unsigned short i = 0; i < PyTuple_Size( pparts ); ++i )
params_.kinematics.minimum_final_state.emplace_back( (PDG)get<int>( PyTuple_GetItem( pparts, i ) ) );
PyObject* pcuts = getElement( kin, "cuts" ); // borrowed
if ( pcuts )
parseParticlesCuts( pcuts );
// for LPAIR/collinear matrix elements
fillLimits( kin, "q2", params_.kinematics.cuts.initial.q2 );
// for the kT factorised matrix elements
fillLimits( kin, "qt", params_.kinematics.cuts.initial.qt );
fillLimits( kin, "phiqt", params_.kinematics.cuts.initial.phi_qt );
fillLimits( kin, "ptdiff", params_.kinematics.cuts.central.pt_diff );
fillLimits( kin, "phiptdiff", params_.kinematics.cuts.central.phi_pt_diff );
fillLimits( kin, "rapiditydiff", params_.kinematics.cuts.central.rapidity_diff );
// generic phase space limits
fillLimits( kin, "rapidity", params_.kinematics.cuts.central.rapidity_single );
fillLimits( kin, "eta", params_.kinematics.cuts.central.eta_single );
fillLimits( kin, "pt", params_.kinematics.cuts.central.pt_single );
fillLimits( kin, "ptsum", params_.kinematics.cuts.central.pt_sum );
fillLimits( kin, "invmass", params_.kinematics.cuts.central.mass_sum );
fillLimits( kin, "mx", params_.kinematics.cuts.remnants.mass_single );
}
void
PythonHandler::parseParticlesCuts( PyObject* cuts )
{
if ( !PyDict_Check( cuts ) )
throwPythonError( "Particle cuts object should be a dictionary!" );
PyObject* pkey = nullptr, *pvalue = nullptr;
Py_ssize_t pos = 0;
while ( PyDict_Next( cuts, &pos, &pkey, &pvalue ) ) {
const PDG pdg = (PDG)get<int>( pkey );
fillLimits( pvalue, "pt", params_.kinematics.cuts.central_particles[pdg].pt_single );
fillLimits( pvalue, "energy", params_.kinematics.cuts.central_particles[pdg].energy_single );
fillLimits( pvalue, "eta", params_.kinematics.cuts.central_particles[pdg].eta_single );
fillLimits( pvalue, "rapidity", params_.kinematics.cuts.central_particles[pdg].rapidity_single );
}
}
void
PythonHandler::parseLogging( PyObject* log )
{
fillParameter( log, "level", (int&)Logger::get().level );
std::vector<std::string> enabled_modules;
fillParameter( log, "enabledModules", enabled_modules );
for ( const auto& mod : enabled_modules )
Logger::get().addExceptionRule( mod );
}
void
PythonHandler::parseIntegrator( PyObject* integr )
{
if ( !PyDict_Check( integr ) )
throwPythonError( "Integrator object should be a dictionary!" );
PyObject* palgo = getElement( integr, MODULE_NAME ); // borrowed
if ( !palgo )
throwPythonError( "Failed to retrieve the integration algorithm name!" );
std::string algo = get<std::string>( palgo );
if ( algo == "plain" )
params_.integrator.type = Integrator::Type::plain;
else if ( algo == "Vegas" ) {
params_.integrator.type = Integrator::Type::Vegas;
fillParameter( integr, "alpha", (double&)params_.integrator.vegas.alpha );
fillParameter( integr, "iterations", params_.integrator.vegas.iterations );
fillParameter( integr, "mode", (int&)params_.integrator.vegas.mode );
fillParameter( integr, "verbosity", (int&)params_.integrator.vegas.verbose );
std::string vegas_logging_output = "cerr";
fillParameter( integr, "loggingOutput", vegas_logging_output );
if ( vegas_logging_output == "cerr" )
// redirect all debugging information to the error stream
params_.integrator.vegas.ostream = stderr;
else if ( vegas_logging_output == "cout" )
// redirect all debugging information to the standard stream
params_.integrator.vegas.ostream = stdout;
else
params_.integrator.vegas.ostream = fopen( vegas_logging_output.c_str(), "w" );
}
else if ( algo == "MISER" ) {
params_.integrator.type = Integrator::Type::MISER;
fillParameter( integr, "estimateFraction", (double&)params_.integrator.miser.estimate_frac );
fillParameter( integr, "minCalls", params_.integrator.miser.min_calls );
fillParameter( integr, "minCallsPerBisection", params_.integrator.miser.min_calls_per_bisection );
fillParameter( integr, "alpha", (double&)params_.integrator.miser.alpha );
fillParameter( integr, "dither", (double&)params_.integrator.miser.dither );
}
else
throwPythonError( Form( "Invalid integration algorithm: %s", algo.c_str() ) );
fillParameter( integr, "numFunctionCalls", params_.integrator.ncvg );
fillParameter( integr, "seed", (unsigned long&)params_.integrator.rng_seed );
unsigned int rng_engine;
fillParameter( integr, "rngEngine", rng_engine );
switch ( rng_engine ) {
case 0: default: params_.integrator.rng_engine = (gsl_rng_type*)gsl_rng_mt19937; break;
case 1: params_.integrator.rng_engine = (gsl_rng_type*)gsl_rng_taus2; break;
case 2: params_.integrator.rng_engine = (gsl_rng_type*)gsl_rng_gfsr4; break;
case 3: params_.integrator.rng_engine = (gsl_rng_type*)gsl_rng_ranlxs0; break;
}
fillParameter( integr, "chiSqCut", params_.integrator.vegas_chisq_cut );
}
void
PythonHandler::parseGenerator( PyObject* gen )
{
if ( !PyDict_Check( gen ) )
throwPythonError( "Generation information object should be a dictionary!" );
params_.generation.enabled = true;
fillParameter( gen, "treat", params_.generation.treat );
fillParameter( gen, "numEvents", params_.generation.maxgen );
fillParameter( gen, "printEvery", params_.generation.gen_print_every );
fillParameter( gen, "numThreads", params_.generation.num_threads );
fillParameter( gen, "numPoints", params_.generation.num_points );
}
void
PythonHandler::parseTamingFunctions( PyObject* tf )
{
if ( !PyList_Check( tf ) )
throwPythonError( "Taming functions list should be a list!" );
for ( Py_ssize_t i = 0; i < PyList_Size( tf ); ++i ) {
PyObject* pit = PyList_GetItem( tf, i ); // borrowed
if ( !pit )
continue;
if ( !PyDict_Check( pit ) )
throwPythonError( Form( "Item %d has invalid type %s", i, pit->ob_type->tp_name ) );
PyObject* pvar = getElement( pit, "variable" ), *pexpr = getElement( pit, "expression" ); // borrowed
params_.taming_functions->add( get<std::string>( pvar ).c_str(), get<std::string>( pexpr ).c_str() );
}
}
void
PythonHandler::parseHadroniser( PyObject* hadr )
{
if ( !PyDict_Check( hadr ) )
throwPythonError( "Hadroniser object should be a dictionary!" );
PyObject* pname = getElement( hadr, MODULE_NAME ); // borrowed
if ( !pname )
throwPythonError( "Hadroniser name is required!" );
std::string hadr_name = get<std::string>( pname );
//--- list of module-specific parameters
ParametersList mod_params;
fillParameter( hadr, "moduleParameters", mod_params );
if ( hadr_name == "pythia8" )
- params_.setHadroniser( new Hadroniser::Pythia8Hadroniser( params_, mod_params ) );
+ params_.setHadroniser( new hadroniser::Pythia8Hadroniser( params_, mod_params ) );
else
throwPythonError( Form( "Unrecognised hadronisation algorithm: \"%s\"!", hadr_name.c_str() ) );
auto h = params_.hadroniser();
{ //--- before calling the init() method
std::vector<std::string> config;
fillParameter( hadr, "preConfiguration", config );
h->readStrings( config );
}
h->init();
{ //--- after init() has been called
std::vector<std::string> config;
fillParameter( hadr, "processConfiguration", config );
for ( const auto& block : config ) {
std::vector<std::string> config_blk;
fillParameter( hadr, block.c_str(), config_blk );
h->readStrings( config_blk );
}
}
}
}
}
#endif
diff --git a/CepGen/Core/Integrand.cpp b/CepGen/Core/Integrand.cpp
index 678cea4..2e4f9ad 100644
--- a/CepGen/Core/Integrand.cpp
+++ b/CepGen/Core/Integrand.cpp
@@ -1,202 +1,201 @@
#include "CepGen/Core/Timer.h"
#include "CepGen/Core/Exception.h"
#include "CepGen/Core/TamingFunction.h"
#include "CepGen/Event/Event.h"
#include "CepGen/Event/Particle.h"
#include "CepGen/Physics/Kinematics.h"
#include "CepGen/Physics/PDG.h"
#include "CepGen/Processes/GenericProcess.h"
#include "CepGen/Hadronisers/GenericHadroniser.h"
#include "CepGen/Parameters.h"
#include <sstream>
#include <fstream>
namespace CepGen
{
namespace Integrand
{
Logger::Level log_level;
Timer tmr;
double
eval( double* x, size_t ndim, void* params )
{
log_level = Logger::get().level;
std::shared_ptr<Event> ev;
Parameters* p = static_cast<Parameters*>( params );
if ( !p )
throw CG_FATAL( "Integrand" ) << "Failed to retrieve the run parameters!";
- Process::GenericProcess* proc = p->process();
+ process::GenericProcess* proc = p->process();
if ( !proc )
throw CG_FATAL( "Integrand" ) << "Failed to retrieve the process!";
//=============================================================================================
// start the timer
//=============================================================================================
tmr.reset();
//=============================================================================================
// prepare the event content prior to the process generation
//=============================================================================================
if ( proc->hasEvent() ) {
ev = proc->event();
if ( proc->first_run ) {
CG_DEBUG( "Integrand" )
<< "Computation launched for " << p->processName() << " process "
<< "0x" << std::hex << p->process() << std::dec << ".\n\t"
<< "Process mode considered: " << p->kinematics.mode << "\n\t"
<< " pz(p1) = " << p->kinematics.incoming_beams.first.pz << "\n\t"
<< " pz(p2) = " << p->kinematics.incoming_beams.second.pz << "\n\t"
<< " structure functions: " << p->kinematics.structure_functions;
p->clearRunStatistics();
proc->first_run = false;
} // passed the first-run preparation
proc->clearEvent();
} // event is not empty
//=============================================================================================
// specify the phase space point to probe
//=============================================================================================
proc->setPoint( ndim, x );
//=============================================================================================
// from this step on, the phase space point is supposed to be set
//=============================================================================================
p->process()->beforeComputeWeight();
double integrand = p->process()->computeWeight();
//=============================================================================================
// invalidate any unphysical behaviour
//=============================================================================================
if ( integrand <= 0. )
return 0.;
//=============================================================================================
// speed up the integration process if no event needs to be generated
//=============================================================================================
if ( !p->storage()
&& !p->taming_functions
&& !p->hadroniser()
&& p->kinematics.cuts.central_particles.size() == 0 )
return integrand;
//=============================================================================================
// fill in the process' Event object
//=============================================================================================
p->process()->fillKinematics();
//=============================================================================================
// once the kinematics variables have been populated, can apply the collection of taming functions
//=============================================================================================
if ( p->taming_functions ) {
if ( p->taming_functions->has( "m_central" )
|| p->taming_functions->has( "pt_central" ) ) {
// build the kinematics of the central system
Particle::Momentum central_system;
for ( const auto& part : ev->getByRole( Particle::CentralSystem ) )
central_system += part.momentum();
// tame the cross-section by the reweighting function
if ( p->taming_functions->has( "m_central" ) )
integrand *= p->taming_functions->eval( "m_central", central_system.mass() );
if ( p->taming_functions->has( "pt_central" ) )
integrand *= p->taming_functions->eval( "pt_central", central_system.pt() );
}
if ( p->taming_functions->has( "q2" ) ) {
integrand *= p->taming_functions->eval( "q2", -ev->getOneByRole( Particle::Parton1 ).momentum().mass() );
integrand *= p->taming_functions->eval( "q2", -ev->getOneByRole( Particle::Parton2 ).momentum().mass() );
}
}
if ( integrand <= 0. )
return 0.;
//=============================================================================================
// set the CepGen part of the event generation
//=============================================================================================
if ( p->storage() )
ev->time_generation = tmr.elapsed();
//=============================================================================================
// event hadronisation and resonances decay
//=============================================================================================
if ( p->hadroniser() ) {
double br = -1.;
if ( !p->hadroniser()->run( *ev, br, p->storage() ) || br == 0. )
return 0.;
integrand *= br; // branching fraction for all decays
}
//=============================================================================================
// apply cuts on final state system (after hadronisation!)
// (watch out your cuts, as this might be extremely time-consuming...)
//=============================================================================================
if ( p->kinematics.cuts.central_particles.size() > 0 ) {
for ( const auto& part : ev->getByRole( Particle::CentralSystem ) ) {
// retrieve all cuts associated to this final state particle
if ( p->kinematics.cuts.central_particles.count( part.pdgId() ) == 0 )
continue;
const auto& cuts_pdgid = p->kinematics.cuts.central_particles.at( part.pdgId() );
// apply these cuts on the given particle
if ( !cuts_pdgid.pt_single.passes( part.momentum().pt() ) )
return 0.;
if ( !cuts_pdgid.energy_single.passes( part.momentum().energy() ) )
return 0.;
if ( !cuts_pdgid.eta_single.passes( part.momentum().eta() ) )
return 0.;
if ( !cuts_pdgid.rapidity_single.passes( part.momentum().rapidity() ) )
return 0.;
}
}
//=============================================================================================
// store the last event in the parameters for its usage by the end user
//=============================================================================================
if ( p->storage() ) {
p->process()->last_event = ev;
p->process()->last_event->time_total = tmr.elapsed();
CG_DEBUG( "Integrand" )
<< "[process 0x" << std::hex << p->process() << std::dec << "] "
<< "Individual time (gen+hadr+cuts): " << p->process()->last_event->time_total*1.e3 << " ms";
}
//=============================================================================================
// a bit of useful debugging
//=============================================================================================
if ( CG_EXCEPT_MATCH( "Integrand", debugInsideLoop ) ) {
std::ostringstream oss;
for ( unsigned short i = 0; i < ndim; ++i )
oss << Form( "%10.8f ", x[i] );
CG_DEBUG( "Integrand" )
<< "f value for dim-" << ndim << " point ( " << oss.str() << "): "
<< integrand;
}
return integrand;
}
}
}
-
diff --git a/CepGen/Core/Parameters.cpp b/CepGen/Core/Parameters.cpp
index a2fc6cf..b1bca39 100644
--- a/CepGen/Core/Parameters.cpp
+++ b/CepGen/Core/Parameters.cpp
@@ -1,257 +1,257 @@
#include "CepGen/Parameters.h"
#include "CepGen/Core/ParametersList.h"
#include "CepGen/Core/Exception.h"
#include "CepGen/Core/TamingFunction.h"
#include "CepGen/Physics/PDG.h"
#include "CepGen/Processes/GenericProcess.h"
#include "CepGen/Hadronisers/GenericHadroniser.h"
#include "CepGen/StructureFunctions/StructureFunctions.h"
#include <iomanip>
namespace CepGen
{
Parameters::Parameters() :
general( new ParametersList ),
taming_functions( new TamingFunctionsCollection ),
store_( false ), total_gen_time_( 0. ), num_gen_events_( 0 )
{}
Parameters::Parameters( Parameters& param ) :
general( param.general ),
kinematics( param.kinematics ), integrator( param.integrator ), generation( param.generation ),
taming_functions( param.taming_functions ),
process_( std::move( param.process_ ) ),
hadroniser_( std::move( param.hadroniser_ ) ),
store_( false ), total_gen_time_( param.total_gen_time_ ), num_gen_events_( param.num_gen_events_ )
{}
Parameters::Parameters( const Parameters& param ) :
general( param.general ),
kinematics( param.kinematics ), integrator( param.integrator ), generation( param.generation ),
taming_functions( param.taming_functions ),
store_( false ), total_gen_time_( param.total_gen_time_ ), num_gen_events_( param.num_gen_events_ )
{}
Parameters::~Parameters() // required for unique_ptr initialisation!
{}
void
Parameters::setThetaRange( float thetamin, float thetamax )
{
kinematics.cuts.central.eta_single = {
Particle::thetaToEta( thetamax ),
Particle::thetaToEta( thetamin )
};
CG_DEBUG( "Parameters" )
<< "eta in range: " << kinematics.cuts.central.eta_single
<< " => theta(min) = " << thetamin << ", theta(max) = " << thetamax << ".";
}
void
Parameters::clearRunStatistics()
{
total_gen_time_ = 0.;
num_gen_events_ = 0;
}
void
Parameters::addGenerationTime( double gen_time )
{
total_gen_time_ += gen_time;
num_gen_events_++;
}
- Process::GenericProcess*
+ process::GenericProcess*
Parameters::process()
{
return process_.get();
}
std::string
Parameters::processName() const
{
if ( !process_ )
return "no process";
return process_->name();
}
void
- Parameters::setProcess( std::unique_ptr<Process::GenericProcess> proc )
+ Parameters::setProcess( std::unique_ptr<process::GenericProcess> proc )
{
process_ = std::move( proc );
}
void
- Parameters::setProcess( Process::GenericProcess* proc )
+ Parameters::setProcess( process::GenericProcess* proc )
{
if ( !proc )
throw CG_FATAL( "Parameters" )
<< "Trying to clone an invalid process!";
process_.reset( proc );
}
- Hadroniser::GenericHadroniser*
+ hadroniser::GenericHadroniser*
Parameters::hadroniser()
{
return hadroniser_.get();
}
std::string
Parameters::hadroniserName() const
{
if ( !hadroniser_ )
return "";
return hadroniser_->name();
}
void
- Parameters::setHadroniser( Hadroniser::GenericHadroniser* hadr )
+ Parameters::setHadroniser( hadroniser::GenericHadroniser* hadr )
{
hadroniser_.reset( hadr );
}
std::ostream&
operator<<( std::ostream& os, const Parameters* p )
{
const bool pretty = true;
const int wb = 90, wt = 40;
os
<< "Parameters dump" << std::left << "\n\n"
<< std::setfill('_') << std::setw( wb+3 ) << "_/¯¯RUN¯INFORMATION¯¯\\_" << std::setfill( ' ' ) << "\n"
<< std::right << std::setw( wb ) << std::left << std::endl
<< std::setw( wt ) << "Process to generate";
if ( p->process_ ) {
os << ( pretty ? boldify( p->process_->name().c_str() ) : p->process_->name() ) << "\n"
<< std::setw( wt ) << "" << p->process_->description();
}
else
os << ( pretty ? boldify( "no process!" ) : "no process!" );
os
<< "\n"
<< std::setw( wt ) << "Events generation? "
<< ( pretty ? yesno( p->generation.enabled ) : std::to_string( p->generation.enabled ) ) << "\n"
<< std::setw( wt ) << "Number of events to generate"
<< ( pretty ? boldify( p->generation.maxgen ) : std::to_string( p->generation.maxgen ) ) << "\n";
if ( p->generation.num_threads > 1 )
os
<< std::setw( wt ) << "Number of threads" << p->generation.num_threads << "\n";
os
<< std::setw( wt ) << "Number of points to try per bin" << p->generation.num_points << "\n"
<< std::setw( wt ) << "Integrand treatment"
<< ( pretty ? yesno( p->generation.treat ) : std::to_string( p->generation.treat ) ) << "\n"
<< std::setw( wt ) << "Verbosity level " << Logger::get().level << "\n";
if ( p->hadroniser_ ) {
os
<< "\n"
<< std::setfill( '-' ) << std::setw( wb+6 )
<< ( pretty ? boldify( " Hadronisation algorithm " ) : "Hadronisation algorithm" ) << std::setfill( ' ' ) << "\n\n"
<< std::setw( wt ) << "Name"
<< ( pretty ? boldify( p->hadroniser_->name().c_str() ) : p->hadroniser_->name() ) << "\n";
}
os
<< "\n"
<< std::setfill( '-' ) << std::setw( wb+6 )
<< ( pretty ? boldify( " Integration parameters " ) : "Integration parameters" ) << std::setfill( ' ' ) << "\n\n";
std::ostringstream int_algo; int_algo << p->integrator.type;
os
<< std::setw( wt ) << "Integration algorithm"
<< ( pretty ? boldify( int_algo.str().c_str() ) : int_algo.str() ) << "\n"
<< std::setw( wt ) << "Number of function calls" << p->integrator.ncvg << "\n"
<< std::setw( wt ) << "Random number generator seed" << p->integrator.rng_seed << "\n";
if ( p->integrator.rng_engine )
os
<< std::setw( wt ) << "Random number generator engine"
<< p->integrator.rng_engine->name << "\n";
std::ostringstream proc_mode; proc_mode << p->kinematics.mode;
os
<< "\n"
<< std::setfill('_') << std::setw( wb+3 )
<< "_/¯¯EVENTS¯KINEMATICS¯¯\\_" << std::setfill( ' ' ) << "\n\n"
<< std::setw( wt ) << "Incoming particles"
<< p->kinematics.incoming_beams.first << ",\n" << std::setw( wt ) << ""
<< p->kinematics.incoming_beams.second << "\n";
if ( p->kinematics.mode != KinematicsMode::invalid )
os << std::setw( wt ) << "Subprocess mode" << ( pretty ? boldify( proc_mode.str().c_str() ) : proc_mode.str() ) << "\n";
if ( p->kinematics.mode != KinematicsMode::ElasticElastic )
os << std::setw( wt ) << "Structure functions" << *p->kinematics.structure_functions << "\n";
os
<< "\n"
<< std::setfill( '-' ) << std::setw( wb+6 ) << ( pretty ? boldify( " Incoming partons " ) : "Incoming partons" ) << std::setfill( ' ' ) << "\n\n";
for ( const auto& lim : p->kinematics.cuts.initial.list() ) // map(particles class, limits)
if ( lim.second.valid() )
os << std::setw( wt ) << lim.first << lim.second << "\n";
os
<< "\n"
<< std::setfill( '-' ) << std::setw( wb+6 ) << ( pretty ? boldify( " Outgoing central system " ) : "Outgoing central system" ) << std::setfill( ' ' ) << "\n\n";
for ( const auto& lim : p->kinematics.cuts.central.list() )
if ( lim.second.valid() )
os << std::setw( wt ) << lim.first << lim.second << "\n";
if ( p->kinematics.cuts.central_particles.size() > 0 ) {
os << std::setw( wt ) << ( pretty ? boldify( ">>> per-particle cuts:" ) : ">>> per-particle cuts:" ) << "\n";
for ( const auto& part_per_lim : p->kinematics.cuts.central_particles ) {
os << " * all single " << std::setw( wt-3 ) << part_per_lim.first << "\n";
for ( const auto& lim : part_per_lim.second.list() )
if ( lim.second.valid() )
os << " - " << std::setw( wt-5 ) << lim.first << lim.second << "\n";
}
}
os << "\n";
os << std::setfill( '-' ) << std::setw( wb+6 ) << ( pretty ? boldify( " Proton / remnants " ) : "Proton / remnants" ) << std::setfill( ' ' ) << "\n\n";
for ( const auto& lim : p->kinematics.cuts.remnants.list() )
os << std::setw( wt ) << lim.first << lim.second << "\n";
return os;
}
std::ostream&
operator<<( std::ostream& os, const Parameters& p )
{
return os << &p;
}
//-----------------------------------------------------------------------------------------------
Parameters::Integration::Integration() :
type( Integrator::Type::Vegas ), ncvg( 500000 ),
rng_seed( 0 ), rng_engine( (gsl_rng_type*)gsl_rng_mt19937 ),
vegas_chisq_cut( 1.5 ),
result( -1. ), err_result( -1. )
{
const size_t ndof = 10;
{
std::shared_ptr<gsl_monte_vegas_state> tmp_state( gsl_monte_vegas_alloc( ndof ), gsl_monte_vegas_free );
gsl_monte_vegas_params_get( tmp_state.get(), &vegas );
}
{
std::shared_ptr<gsl_monte_miser_state> tmp_state( gsl_monte_miser_alloc( ndof ), gsl_monte_miser_free );
gsl_monte_miser_params_get( tmp_state.get(), &miser );
}
}
Parameters::Integration::Integration( const Integration& rhs ) :
type( rhs.type ), ncvg( rhs.ncvg ),
rng_seed( rhs.rng_seed ), rng_engine( rhs.rng_engine ),
vegas( rhs.vegas ), vegas_chisq_cut( rhs.vegas_chisq_cut ),
miser( rhs.miser ),
result( -1. ), err_result( -1. )
{}
Parameters::Integration::~Integration()
{
//if ( vegas.ostream && vegas.ostream != stdout && vegas.ostream != stderr )
// fclose( vegas.ostream );
}
//-----------------------------------------------------------------------------------------------
Parameters::Generation::Generation() :
enabled( false ), maxgen( 0 ),
symmetrise( false ), treat( false ), ngen( 0 ), gen_print_every( 10000 ),
num_threads( 2 ), num_points( 100 )
{}
}
diff --git a/CepGen/Hadronisers/GenericHadroniser.cpp b/CepGen/Hadronisers/GenericHadroniser.cpp
index b2049da..8b297ec 100644
--- a/CepGen/Hadronisers/GenericHadroniser.cpp
+++ b/CepGen/Hadronisers/GenericHadroniser.cpp
@@ -1,54 +1,53 @@
#include "CepGen/Hadronisers/GenericHadroniser.h"
#include "CepGen/Core/Exception.h"
#include "CepGen/Core/ParametersList.h"
namespace CepGen
{
- namespace Hadroniser
+ namespace hadroniser
{
GenericHadroniser::GenericHadroniser( const char* name, const ParametersList& plist ) :
name_( name ),
seed_ ( plist.get<int>( "seed", -1ll ) ),
max_trials_( plist.get<int>( "maxTrials", 1 ) )
{
CG_DEBUG( "Hadroniser:init" )
<< "\"" << name_ << "\"-type hadroniser built with:\n\t"
<< "* seed = " << seed_ << "\n\t"
<< "* maximum trials: " << max_trials_;
}
void
GenericHadroniser::readStrings( const std::vector<std::string>& params ) {
if ( params.empty() )
return;
std::ostringstream os;
for ( const auto& p : params ) {
readString( p );
os << "\n\t '" << p << "'";
}
CG_DEBUG( "Hadroniser:configure" )
<< "Feeding \"" << name_ << "\" hadroniser with:"
<< os.str();
}
std::string
GenericHadroniser::name() const
{
return name_;
}
}
std::ostream&
- operator<<( std::ostream& os, const Hadroniser::GenericHadroniser& hadr )
+ operator<<( std::ostream& os, const hadroniser::GenericHadroniser& hadr )
{
return os << hadr.name().c_str();
}
std::ostream&
- operator<<( std::ostream& os, const Hadroniser::GenericHadroniser* hadr )
+ operator<<( std::ostream& os, const hadroniser::GenericHadroniser* hadr )
{
return os << hadr->name().c_str();
}
}
-
diff --git a/CepGen/Hadronisers/GenericHadroniser.h b/CepGen/Hadronisers/GenericHadroniser.h
index b7b7440..082fe9c 100644
--- a/CepGen/Hadronisers/GenericHadroniser.h
+++ b/CepGen/Hadronisers/GenericHadroniser.h
@@ -1,70 +1,70 @@
#ifndef CepGen_Hadronisers_GenericHadroniser_h
#define CepGen_Hadronisers_GenericHadroniser_h
#include <vector>
#include <memory>
#include <iostream>
namespace CepGen
{
class Event;
class Particle;
class Parameters;
class ParametersList;
/// Location for all hadronisers to be run downstream to the events generation
- namespace Hadroniser
+ namespace hadroniser
{
/**
* \brief Class template to define any hadroniser as a general object with defined methods
* \author Laurent Forthomme <laurent.forthomme@cern.ch>
* \date January 2014
*/
class GenericHadroniser
{
public:
/// Write out all hadroniser attributes in output stream
friend std::ostream& operator<<( std::ostream& os, const GenericHadroniser& hadr );
/// Write out all hadroniser attributes in output stream
friend std::ostream& operator<<( std::ostream& os, const GenericHadroniser* hadr );
/// Default constructor for an undefined hadroniser
explicit GenericHadroniser( const char* name, const ParametersList& );
virtual ~GenericHadroniser() {}
/// Parse a configuration string
virtual void readString( const char* ) {}
/// Parse a configuration string
virtual void readString( const std::string& param ) { readString( param.c_str() ); }
/// Parse a list of configuration strings
virtual void readStrings( const std::vector<std::string>& params );
/// Initialise the event hadroniser before its running
virtual void init() = 0;
/** \brief Hadronise a full event
* \param[inout] ev Event to hadronise
* \param[inout] weight Event weight after hadronisation
* \param[in] full Perform the full state hadronisation (incl. remnants fragmentation)
* \return Boolean stating whether or not the hadronisation occured successfully
*/
virtual bool run( Event& ev, double& weight, bool full ) = 0;
/// Specify the process cross section, in pb
virtual void setCrossSection( double xsec, double xsec_err ) {}
/// \brief Specify a random numbers generator seed for the hadroniser
/// \param[in] seed A RNG seed
void setSeed( long long seed ) { seed_ = seed; }
/// Return a human-readable name for this hadroniser
std::string name() const;
protected:
/// Name of the hadroniser
std::string name_;
/// Random numbers generator seed for the hadroniser
long long seed_;
/// Maximal number of trials for the hadronisation of the proton(s) remnants
unsigned short max_trials_;
};
}
}
#endif
diff --git a/CepGen/Hadronisers/Pythia8Hadroniser.cpp b/CepGen/Hadronisers/Pythia8Hadroniser.cpp
index 85964dd..ca6afb6 100644
--- a/CepGen/Hadronisers/Pythia8Hadroniser.cpp
+++ b/CepGen/Hadronisers/Pythia8Hadroniser.cpp
@@ -1,464 +1,464 @@
#include "CepGen/Hadronisers/Pythia8Hadroniser.h"
#include "CepGen/Parameters.h"
#include "CepGen/Physics/Kinematics.h"
#include "CepGen/Physics/Constants.h"
#include "CepGen/Physics/PDG.h"
#include "CepGen/Core/ParametersList.h"
#include "CepGen/Core/Exception.h"
#include "CepGen/Core/utils.h"
#include "CepGen/Event/Event.h"
#include "CepGen/Event/Particle.h"
#include "CepGen/Version.h"
namespace CepGen
{
- namespace Hadroniser
+ namespace hadroniser
{
#ifdef PYTHIA8
Pythia8::Vec4
momToVec4( const Particle::Momentum& mom )
{
return Pythia8::Vec4( mom.px(), mom.py(), mom.pz(), mom.energy() );
}
#endif
Pythia8Hadroniser::Pythia8Hadroniser( const Parameters& params, const ParametersList& plist ) :
GenericHadroniser( "pythia8", plist ),
#ifdef PYTHIA8
pythia_( new Pythia8::Pythia ), lhaevt_( new LHAEvent( ¶ms ) ),
#endif
full_evt_( false ), offset_( 0 ), first_evt_( true ), params_( ¶ms )
{
#ifdef PYTHIA8
pythia_->setLHAupPtr( (Pythia8::LHAup*)lhaevt_.get() );
pythia_->settings.parm( "Beams:idA", (short)params.kinematics.incoming_beams.first.pdg );
pythia_->settings.parm( "Beams:idB", (short)params.kinematics.incoming_beams.second.pdg );
// specify we will be using a LHA input
pythia_->settings.mode( "Beams:frameType", 5 );
pythia_->settings.parm( "Beams:eCM", params.kinematics.sqrtS() );
#endif
for ( const auto& pdgid : params.kinematics.minimum_final_state )
min_ids_.emplace_back( (unsigned short)pdgid );
}
Pythia8Hadroniser::~Pythia8Hadroniser()
{
#ifdef PYTHIA8
pythia_->settings.writeFile( "last_pythia_config.cmd", false );
#endif
}
void
Pythia8Hadroniser::readString( const char* param )
{
#ifdef PYTHIA8
if ( !pythia_->readString( param ) )
throw CG_FATAL( "Pythia8Hadroniser" ) << "The Pythia8 core failed to parse the following setting:\n\t" << param;
#endif
}
void
Pythia8Hadroniser::init()
{
#ifdef PYTHIA8
if ( pythia_->settings.flag( "ProcessLevel:all" ) != full_evt_ )
pythia_->settings.flag( "ProcessLevel:all", full_evt_ );
if ( seed_ == -1ll )
pythia_->settings.flag( "Random:setSeed", false );
else {
pythia_->settings.flag( "Random:setSeed", true );
pythia_->settings.mode( "Random:seed", seed_ );
}
switch ( params_->kinematics.mode ) {
case KinematicsMode::ElasticElastic: {
pythia_->settings.mode( "BeamRemnants:unresolvedHadron", 3 );
} break;
case KinematicsMode::InelasticElastic: {
pythia_->settings.mode( "BeamRemnants:unresolvedHadron", 2 );
} break;
case KinematicsMode::ElasticInelastic: {
pythia_->settings.mode( "BeamRemnants:unresolvedHadron", 1 );
} break;
case KinematicsMode::InelasticInelastic: default: {
pythia_->settings.mode( "BeamRemnants:unresolvedHadron", 0 );
} break;
}
if ( !pythia_->init() )
throw CG_FATAL( "Pythia8Hadroniser" )
<< "Failed to initialise the Pythia8 core!\n\t"
<< "See the message above for more details.";
#else
throw CG_FATAL( "Pythia8Hadroniser" )
<< "Pythia8 is not linked to this instance!";
#endif
}
void
Pythia8Hadroniser::setCrossSection( double xsec, double xsec_err )
{
#ifdef PYTHIA8
lhaevt_->setCrossSection( 0, xsec, xsec_err );
#endif
}
bool
Pythia8Hadroniser::run( Event& ev, double& weight, bool full )
{
//--- initialise the event weight before running any decay algorithm
weight = 1.;
#ifdef PYTHIA8
if ( !full && !pythia_->settings.flag( "ProcessLevel:resonanceDecays" ) )
return true;
//--- switch full <-> partial event
if ( full != full_evt_ ) {
full_evt_ = full;
init();
}
//===========================================================================================
// convert our event into a custom LHA format
//===========================================================================================
lhaevt_->feedEvent( ev, full, params_->kinematics.mode );
//if ( full ) lhaevt_->listEvent();
//===========================================================================================
// launch the hadronisation / resonances decays, and update the event accordingly
//===========================================================================================
ev.num_hadronisation_trials = 0;
while ( true ) {
if ( ev.num_hadronisation_trials++ > max_trials_ )
return false;
//--- run the hadronisation/fragmentation algorithm
if ( pythia_->next() ) {
//--- hadronisation successful
if ( first_evt_ && full ) {
offset_ = 0;
for ( unsigned short i = 1; i < pythia_->event.size(); ++i )
if ( pythia_->event[i].status() == -12 ) // skip the incoming particles
offset_++;
first_evt_ = false;
}
break;
}
}
//===========================================================================================
// update the event content with Pythia's output
//===========================================================================================
updateEvent( ev, weight, full );
return true;
#else
throw CG_FATAL( "Pythia8Hadroniser" ) << "Pythia8 is not linked to this instance!";
#endif
}
#ifdef PYTHIA8
Particle&
Pythia8Hadroniser::addParticle( Event& ev, const Pythia8::Particle& py_part, const Pythia8::Vec4& mom, unsigned short role ) const
{
Particle& op = ev.addParticle( (Particle::Role)role );
op.setPdgId( static_cast<PDG>( abs( py_part.id() ) ), py_part.charge() );
op.setStatus( py_part.isFinal()
? Particle::Status::FinalState
: Particle::Status::Propagator );
op.setMomentum( Particle::Momentum( mom.px(), mom.py(), mom.pz(), mom.e() ) );
op.setMass( mom.mCalc() );
lhaevt_->addCorresp( py_part.index()-offset_, op.id() );
return op;
}
void
Pythia8Hadroniser::updateEvent( Event& ev, double& weight, bool full ) const
{
for ( unsigned short i = 1+offset_; i < pythia_->event.size(); ++i ) {
const Pythia8::Particle& p = pythia_->event[i];
const unsigned short cg_id = lhaevt_->cepgenId( i-offset_ );
if ( cg_id != LHAEvent::invalid_id ) {
//----- particle already in the event
Particle& cg_part = ev[cg_id];
//--- fragmentation result
if ( cg_part.role() == Particle::OutgoingBeam1
|| cg_part.role() == Particle::OutgoingBeam2 ) {
cg_part.setStatus( Particle::Status::Fragmented );
continue;
}
//--- particle is not what we expect
if ( p.idAbs() != abs( cg_part.integerPdgId() ) ) {
CG_INFO( "Pythia8Hadroniser:update" ) << "LHAEVT event content:";
lhaevt_->listEvent();
CG_INFO( "Pythia8Hadroniser:update" ) << "Pythia event content:";
pythia_->event.list();
CG_INFO( "Pythia8Hadroniser:update" ) << "CepGen event content:";
ev.dump();
CG_INFO( "Pythia8Hadroniser:update" ) << "Correspondence:";
lhaevt_->dumpCorresp();
throw CG_FATAL( "Pythia8Hadroniser:update" )
<< "Event list corruption detected for (Pythia/CepGen) particle " << i << "/" << cg_id << ":\n\t"
<< "should be " << abs( p.id() ) << ", "
<< "got " << cg_part.integerPdgId() << "!";
}
//--- resonance decayed; apply branching ratio for this decay
if ( p.particleDataEntry().sizeChannels() > 0 ) {
weight *= p.particleDataEntry().pickChannel().bRatio();
cg_part.setStatus( Particle::Status::Resonance );
}
}
else {
//----- new particle to be added
const unsigned short role = findRole( ev, p );
switch ( (Particle::Role)role ) {
default: break;
case Particle::OutgoingBeam1: {
ev.getByRole( Particle::OutgoingBeam1 )[0].setStatus( Particle::Status::Fragmented );
if ( abs( p.status() ) != 61 )
break;
} // no break!
case Particle::OutgoingBeam2: {
ev.getByRole( Particle::OutgoingBeam2 )[0].setStatus( Particle::Status::Fragmented );
if ( abs( p.status() ) != 61 )
break;
} // no break!
}
// found the role ; now we can add the particle
Particle& cg_part = addParticle( ev, p, p.p(), role );
for ( const auto& moth_id : p.motherList() ) {
if ( moth_id <= offset_ )
continue;
const unsigned short moth_cg_id = lhaevt_->cepgenId( moth_id-offset_ );
if ( moth_cg_id != LHAEvent::invalid_id )
cg_part.addMother( ev[moth_cg_id] );
else
cg_part.addMother( addParticle( ev, pythia_->event[moth_id], p.p(), role ) );
if ( !p.isFinal() ) {
if ( p.isResonance() || p.daughterList().size() > 0 )
cg_part.setStatus( Particle::Status::Resonance );
else
cg_part.setStatus( Particle::Status::Undefined );
}
}
}
}
}
unsigned short
Pythia8Hadroniser::findRole( const Event& ev, const Pythia8::Particle& p ) const
{
for ( const auto& par_id : p.motherList() ) {
if ( par_id == 1 && offset_ > 0 )
return (unsigned short)Particle::OutgoingBeam1;
if ( par_id == 2 && offset_ > 0 )
return (unsigned short)Particle::OutgoingBeam2;
const unsigned short par_cg_id = lhaevt_->cepgenId( par_id-offset_ );
if ( par_cg_id != LHAEvent::invalid_id )
return (unsigned short)ev.at( par_cg_id ).role();
return findRole( ev, pythia_->event[par_id] );
}
return (unsigned short)Particle::UnknownRole;
}
#endif
}
//================================================================================================
// Custom LHA event definition
//================================================================================================
#ifdef PYTHIA8
const double LHAEvent::mp_ = ParticleProperties::mass( PDG::proton );
const double LHAEvent::mp2_ = LHAEvent::mp_*LHAEvent::mp_;
LHAEvent::LHAEvent( const Parameters* params ) :
LHAup( 3 ), params_( params )
{
addProcess( 0, 1., 1., 1.e3 );
if ( params_ ) {
setBeamA( (short)params_->kinematics.incoming_beams.first.pdg, params_->kinematics.incoming_beams.first.pz );
setBeamB( (short)params_->kinematics.incoming_beams.second.pdg, params_->kinematics.incoming_beams.second.pz );
}
}
void
LHAEvent::setCrossSection( int id, double xsec, double xsec_err )
{
setXSec( id, xsec );
setXErr( id, xsec_err );
//listInit();
}
void
LHAEvent::feedEvent( const Event& ev, bool full, const KinematicsMode& mode )
{
const double scale = ev.getOneByRole( Particle::Intermediate ).mass();
setProcess( 0, 1., scale, Constants::alphaEM, Constants::alphaQCD );
const Particle& part1 = ev.getOneByRole( Particle::Parton1 ), &part2 = ev.getOneByRole( Particle::Parton2 );
const Particle& op1 = ev.getOneByRole( Particle::OutgoingBeam1 ), &op2 = ev.getOneByRole( Particle::OutgoingBeam2 );
const double q2_1 = -part1.momentum().mass2(), q2_2 = -part2.momentum().mass2();
const double x1 = q2_1/( q2_1+op1.mass2()-mp2_ ), x2 = q2_2/( q2_2+op2.mass2()-mp2_ );
unsigned short quark1_id = 0, quark2_id = 0;
unsigned short quark1_pdgid = part1.integerPdgId(), quark2_pdgid = part2.integerPdgId();
- const Pythia8::Vec4 mom_part1( Hadroniser::momToVec4( part1.momentum() ) ), mom_part2( Hadroniser::momToVec4( part2.momentum() ) );
+ const Pythia8::Vec4 mom_part1( hadroniser::momToVec4( part1.momentum() ) ), mom_part2( hadroniser::momToVec4( part2.momentum() ) );
if ( !full ) {
//=============================================================================================
// incoming partons
//=============================================================================================
addCorresp( sizePart(), part1.id() );
addParticle( quark1_pdgid, -2, quark1_id, 0, 0, 0, mom_part1.px(), mom_part1.py(), mom_part1.pz(), mom_part1.e(), mom_part1.mCalc(), 0., 0. );
addCorresp( sizePart(), part2.id() );
addParticle( quark2_pdgid, -2, quark2_id, 0, 0, 0, mom_part2.px(), mom_part2.py(), mom_part2.pz(), mom_part2.e(), mom_part2.mCalc(), 0., 0. );
}
else { // full event content (with collinear partons)
const bool inel1 = ( mode == KinematicsMode::InelasticElastic || mode == KinematicsMode::InelasticInelastic );
const bool inel2 = ( mode == KinematicsMode::ElasticInelastic || mode == KinematicsMode::InelasticInelastic );
Pythia8::Vec4 mom_iq1 = mom_part1, mom_iq2 = mom_part2;
unsigned short colour_index = 501, quark1_colour = 0, quark2_colour = 0;
//FIXME select quark flavours accordingly
if ( inel1 ) {
quark1_pdgid = 2;
quark1_colour = colour_index++;
- mom_iq1 = Hadroniser::momToVec4( x1*ev.getOneByRole( Particle::IncomingBeam1 ).momentum() );
+ mom_iq1 = hadroniser::momToVec4( x1*ev.getOneByRole( Particle::IncomingBeam1 ).momentum() );
}
if ( inel2 ) {
quark2_pdgid = 2;
quark2_colour = colour_index++;
- mom_iq2 = Hadroniser::momToVec4( x2*ev.getOneByRole( Particle::IncomingBeam2 ).momentum() );
+ mom_iq2 = hadroniser::momToVec4( x2*ev.getOneByRole( Particle::IncomingBeam2 ).momentum() );
}
//--- flavour / x value of hard-process initiators
setIdX( part1.integerPdgId(), part2.integerPdgId(), x1, x2 );
//===========================================================================================
// incoming valence quarks
//===========================================================================================
quark1_id = sizePart();
addCorresp( quark1_id, op1.id() );
addParticle( quark1_pdgid, -1, 0, 0, quark1_colour, 0, mom_iq1.px(), mom_iq1.py(), mom_iq1.pz(), mom_iq1.e(), mom_iq1.mCalc(), 0., 1. );
quark2_id = sizePart();
addCorresp( quark2_id, op2.id() );
addParticle( quark2_pdgid, -1, 0, 0, quark2_colour, 0, mom_iq2.px(), mom_iq2.py(), mom_iq2.pz(), mom_iq2.e(), mom_iq2.mCalc(), 0., 1. );
//===========================================================================================
// outgoing valence quarks
//===========================================================================================
if ( inel1 ) {
const Pythia8::Vec4 mom_oq1 = mom_iq1-mom_part1;
addParticle( quark1_pdgid, 1, quark1_id, quark2_id, quark1_colour, 0, mom_oq1.px(), mom_oq1.py(), mom_oq1.pz(), mom_oq1.e(), mom_oq1.mCalc(), 0., 1. );
}
if ( inel2 ) {
const Pythia8::Vec4 mom_oq2 = mom_iq2-mom_part2;
addParticle( quark2_pdgid, 1, quark1_id, quark2_id, quark2_colour, 0, mom_oq2.px(), mom_oq2.py(), mom_oq2.pz(), mom_oq2.e(), mom_oq2.mCalc(), 0., 1. );
}
}
//=============================================================================================
// central system
//=============================================================================================
for ( const auto& p : ev.getByRole( Particle::CentralSystem ) ) {
const auto mothers = p.mothers();
unsigned short moth1_id = 1, moth2_id = 2;
if ( !full ) {
moth1_id = moth2_id = 0;
if ( mothers.size() > 0 ) {
const unsigned short moth1_cg_id = *mothers.begin();
moth1_id = pythiaId( moth1_cg_id );
if ( moth1_id == invalid_id ) {
const Particle& moth = ev.at( moth1_cg_id );
if ( moth.mothers().size() > 0 )
moth1_id = pythiaId( *moth.mothers().begin() );
if ( moth.mothers().size() > 1 )
moth2_id = pythiaId( *moth.mothers().rbegin() );
}
if ( mothers.size() > 1 ) {
const unsigned short moth2_cg_id = *mothers.rbegin();
moth2_id = pythiaId( moth2_cg_id );
if ( moth2_id == invalid_id ) {
const Particle& moth = ev.at( moth2_cg_id );
moth.dump();
moth2_id = pythiaId( *moth.mothers().rbegin() );
}
}
}
}
const Pythia8::Vec4 mom_part( p.momentum().px(), p.momentum().py(), p.momentum().pz(), p.momentum().energy() );
addCorresp( sizePart(), p.id() );
addParticle( p.integerPdgId(), 1, moth1_id, moth2_id, 0, 0, mom_part.px(), mom_part.py(), mom_part.pz(), mom_part.e(), mom_part.mCalc(), 0., 0., 0. );
}
setPdf( quark1_pdgid, quark2_pdgid, x1, x2, scale, 0., 0., false );
}
bool
LHAEvent::setInit()
{
return true;
}
bool
LHAEvent::setEvent( int )
{
return true;
}
void
LHAEvent::setProcess( int id, double xsec, double q2_scale, double alpha_qed, double alpha_qcd )
{
LHAup::setProcess( id, xsec, q2_scale, alpha_qed, alpha_qcd );
py_cg_corresp_.clear();
}
unsigned short
LHAEvent::cepgenId( unsigned short py_id ) const
{
for ( const auto& py_cg : py_cg_corresp_ )
if ( py_cg.first == py_id )
return py_cg.second;
return invalid_id;
}
unsigned short
LHAEvent::pythiaId( unsigned short cg_id ) const
{
for ( const auto& py_cg : py_cg_corresp_ )
if ( py_cg.second == cg_id )
return py_cg.first;
return invalid_id;
}
void
LHAEvent::addCorresp( unsigned short py_id, unsigned short cg_id )
{
py_cg_corresp_.emplace_back( py_id, cg_id );
}
void
LHAEvent::dumpCorresp() const
{
std::ostringstream oss;
oss << "List of Pythia <-> CepGen particle ids correspondance";
for ( const auto& py_cg : py_cg_corresp_ )
oss << "\n\t" << py_cg.first << " <-> " << py_cg.second;
CG_INFO( "LHAEvent:dump" ) << oss.str();
}
#endif
}
diff --git a/CepGen/Hadronisers/Pythia8Hadroniser.h b/CepGen/Hadronisers/Pythia8Hadroniser.h
index dfdbcbc..d4dea88 100644
--- a/CepGen/Hadronisers/Pythia8Hadroniser.h
+++ b/CepGen/Hadronisers/Pythia8Hadroniser.h
@@ -1,84 +1,84 @@
#ifndef CepGen_Hadronisers_Pythia8Hadroniser_h
#define CepGen_Hadronisers_Pythia8Hadroniser_h
#include "CepGen/Hadronisers/GenericHadroniser.h"
#ifdef PYTHIA8
#include <Pythia8/Pythia.h>
#include <memory>
#endif
#include <unordered_map>
#include <vector>
namespace CepGen
{
class Particle;
class ParametersList;
enum class KinematicsMode;
#ifdef PYTHIA8
class LHAEvent : public Pythia8::LHAup
{
public:
explicit LHAEvent( const Parameters* );
void feedEvent( const Event& ev, bool full, const KinematicsMode& );
bool setInit() override;
bool setEvent( int ) override;
void setCrossSection( int id, double xsec, double xsec_err );
void setProcess( int id, double xsec, double q2_scale, double alpha_qed, double alpha_qcd );
unsigned short cepgenId( unsigned short py_id ) const;
unsigned short pythiaId( unsigned short cg_id ) const;
void addCorresp( unsigned short py_id, unsigned short cg_id );
void dumpCorresp() const;
static constexpr unsigned short invalid_id = 999;
private:
static const double mp_, mp2_;
std::vector<std::pair<unsigned short, unsigned short> > py_cg_corresp_;
const Parameters* params_;
};
#endif
- namespace Hadroniser
+ namespace hadroniser
{
/**
* Full interface to the Pythia8 hadronisation algorithm. It can be used in a single particle decay mode as well as a full event hadronisation using the string model, as in Jetset.
* \brief Pythia8 hadronisation algorithm
*/
class Pythia8Hadroniser : public GenericHadroniser
{
public:
explicit Pythia8Hadroniser( const Parameters&, const ParametersList& );
~Pythia8Hadroniser();
void readString( const char* param ) override;
void init() override;
bool run( Event& ev, double& weight, bool full ) override;
void setCrossSection( double xsec, double xsec_err ) override;
bool fullEvent() const { return full_evt_; }
void setFullEvent( bool full = true ) { full_evt_ = full; }
private:
static constexpr unsigned short invalid_idx_ = 999;
std::vector<unsigned short> min_ids_;
std::unordered_map<short,short> py_cg_corresp_;
#ifdef PYTHIA8
unsigned short findRole( const Event& ev, const Pythia8::Particle& p ) const;
void updateEvent( Event& ev, double& weight, bool full ) const;
Particle& addParticle( Event& ev, const Pythia8::Particle&, const Pythia8::Vec4& mom, unsigned short ) const;
/// A Pythia8 core to be wrapped
std::unique_ptr<Pythia8::Pythia> pythia_;
std::unique_ptr<LHAEvent> lhaevt_;
#endif
bool full_evt_;
unsigned short offset_;
bool first_evt_;
const Parameters* params_; // not owning
};
}
}
#endif
diff --git a/CepGen/Parameters.h b/CepGen/Parameters.h
index 2e12481..a9c0e4c 100644
--- a/CepGen/Parameters.h
+++ b/CepGen/Parameters.h
@@ -1,147 +1,147 @@
#ifndef CepGen_Parameters_h
#define CepGen_Parameters_h
#include "CepGen/Core/Integrator.h"
#include "CepGen/Physics/Kinematics.h"
#include <memory>
#include <gsl/gsl_monte_vegas.h>
#include <gsl/gsl_monte_miser.h>
namespace CepGen
{
class Event;
class TamingFunctionsCollection;
class ParametersList;
- namespace Process { class GenericProcess; }
- namespace Hadroniser { class GenericHadroniser; }
+ namespace process { class GenericProcess; }
+ namespace hadroniser { class GenericHadroniser; }
/// List of parameters used to start and run the simulation job
class Parameters
{
public:
Parameters();
/// Copy constructor (transfers ownership to the process/hadroniser!)
Parameters( Parameters& );
/// Const copy constructor (all but the process and the hadroniser)
Parameters( const Parameters& );
~Parameters(); // required for unique_ptr initialisation!
/// Set the polar angle range for the produced leptons
/// \param[in] thetamin The minimal value of \f$\theta\f$ for the outgoing leptons
/// \param[in] thetamax The maximal value of \f$\theta\f$ for the outgoing leptons
void setThetaRange( float thetamin, float thetamax );
/// Dump the input parameters in the terminal
friend std::ostream& operator<<( std::ostream& os, const Parameters* );
friend std::ostream& operator<<( std::ostream& os, const Parameters& );
std::shared_ptr<ParametersList> general;
//----- process to compute
/// Process for which the cross-section will be computed and the events will be generated
- Process::GenericProcess* process();
+ process::GenericProcess* process();
/// Name of the process considered
std::string processName() const;
/// Set the process to study
- void setProcess( std::unique_ptr<Process::GenericProcess> proc );
+ void setProcess( std::unique_ptr<process::GenericProcess> proc );
/// Set the process to study
- void setProcess( Process::GenericProcess* proc );
+ void setProcess( process::GenericProcess* proc );
//----- events kinematics
/// Events kinematics for phase space definition
Kinematics kinematics;
//----- VEGAS
/// Collection of integrator parameters
struct Integration
{
Integration();
Integration( const Integration& );
~Integration();
Integrator::Type type;
/// Number of function calls to be computed for each point
unsigned int ncvg; // ??
/// Random number generator seed
long rng_seed;
/// Random number generator engine
gsl_rng_type* rng_engine;
gsl_monte_vegas_params vegas;
double vegas_chisq_cut;
gsl_monte_miser_params miser;
double result, err_result;
};
/// Integrator parameters
Integration integrator;
//----- events generation
/// Collection of events generation parameters
struct Generation
{
Generation();
/// Are we generating events ? (true) or are we only computing the cross-section ? (false)
bool enabled;
/// Maximal number of events to generate in this run
unsigned int maxgen;
/// Do we want the events to be symmetrised with respect to the \f$z\f$-axis ?
bool symmetrise;
/// Is the integrand to be smoothed for events generation?
bool treat;
/// Number of events already generated in this run
unsigned int ngen;
/// Frequency at which the events are displayed to the end-user
unsigned int gen_print_every;
/// Number of threads to perform the integration
unsigned int num_threads;
/// Number of points to "shoot" in each integration bin by the algorithm
unsigned int num_points;
};
/// Events generation parameters
Generation generation;
/// Specify if the generated events are to be stored
void setStorage( bool store ) { store_ = store; }
/// Are the events generated in this run to be stored in the output file ?
bool storage() const { return store_; }
//----- hadronisation algorithm
/// Hadronisation algorithm to use for the proton(s) fragmentation
- Hadroniser::GenericHadroniser* hadroniser();
+ hadroniser::GenericHadroniser* hadroniser();
/// Name of the hadroniser (if applicable)
std::string hadroniserName() const;
/// Set the hadronisation algorithm
- void setHadroniser( Hadroniser::GenericHadroniser* hadr );
+ void setHadroniser( hadroniser::GenericHadroniser* hadr );
//----- taming functions
/// Functionals to be used to account for rescattering corrections (implemented within the process)
std::shared_ptr<TamingFunctionsCollection> taming_functions;
//----- run statistics
/// Reset the total generation time and the number of events generated for this run
void clearRunStatistics();
/// Add a new timing into the total generation time
/// \param[in] gen_time Time to add (in seconds)
void addGenerationTime( double gen_time );
/// Return the total generation time for this run (in seconds)
inline double totalGenerationTime() const { return total_gen_time_; }
/// Total number of events already generated in this run
inline unsigned int numGeneratedEvents() const { return num_gen_events_; }
private:
- std::unique_ptr<Process::GenericProcess> process_;
- std::unique_ptr<Hadroniser::GenericHadroniser> hadroniser_;
+ std::unique_ptr<process::GenericProcess> process_;
+ std::unique_ptr<hadroniser::GenericHadroniser> hadroniser_;
bool store_;
/// Total generation time (in seconds)
double total_gen_time_;
/// Number of events already generated
unsigned int num_gen_events_;
};
}
#endif
diff --git a/CepGen/Processes/FortranKTProcess.cpp b/CepGen/Processes/FortranKTProcess.cpp
index ca67259..9124417 100644
--- a/CepGen/Processes/FortranKTProcess.cpp
+++ b/CepGen/Processes/FortranKTProcess.cpp
@@ -1,167 +1,166 @@
#include "CepGen/Processes/FortranKTProcess.h"
#include "CepGen/Processes/Fortran/KTStructures.h"
#include "CepGen/Core/ParametersList.h"
#include "CepGen/StructureFunctions/StructureFunctions.h"
#include "CepGen/Event/Event.h"
#include "CepGen/Physics/KTFlux.h"
#include "CepGen/Physics/Constants.h"
#include "CepGen/Physics/PDG.h"
extern "C"
{
extern CepGen::KTBlock::Constants constants_;
extern CepGen::KTBlock::Parameters params_;
extern CepGen::KTBlock::KTKinematics ktkin_;
extern CepGen::KTBlock::Cuts kincuts_;
extern CepGen::KTBlock::Event evtkin_;
}
namespace CepGen
{
- namespace Process
+ namespace process
{
FortranKTProcess::FortranKTProcess( const ParametersList& params, const char* name, const char* descr, std::function<void( double& )> func ) :
GenericKTProcess( params, name, descr, { { PDG::photon, PDG::photon } }, { PDG::muon, PDG::muon } ),
pair_( params.get<int>( "pair", 13 ) ),
method_( params.get<int>( "method", 1 ) ),
func_( func )
{
constants_.m_p = GenericProcess::mp_;
constants_.units = Constants::GeV2toBarn;
constants_.pi = M_PI;
constants_.alpha_em = Constants::alphaEM;
}
void
FortranKTProcess::preparePhaseSpace()
{
mom_ip1_ = event_->getOneByRole( Particle::IncomingBeam1 ).momentum();
mom_ip2_ = event_->getOneByRole( Particle::IncomingBeam2 ).momentum();
registerVariable( y1_, Mapping::linear, cuts_.cuts.central.rapidity_single, { -6., 6. }, "First central particle rapidity" );
registerVariable( y2_, Mapping::linear, cuts_.cuts.central.rapidity_single, { -6., 6. }, "Second central particle rapidity" );
registerVariable( pt_diff_, Mapping::linear, cuts_.cuts.central.pt_diff, { 0., 50. }, "Transverse momentum difference between central particles" );
registerVariable( phi_pt_diff_, Mapping::linear, cuts_.cuts.central.phi_pt_diff, { 0., 2.*M_PI }, "Central particles azimuthal angle difference" );
//===========================================================================================
// feed phase space cuts to the common block
//===========================================================================================
cuts_.cuts.central.pt_single.save( (bool&)kincuts_.ipt, kincuts_.pt_min, kincuts_.pt_max );
cuts_.cuts.central.energy_single.save( (bool&)kincuts_.iene, kincuts_.ene_min, kincuts_.ene_max );
cuts_.cuts.central.eta_single.save( (bool&)kincuts_.ieta, kincuts_.eta_min, kincuts_.eta_max );
cuts_.cuts.central.mass_sum.save( (bool&)kincuts_.iinvm, kincuts_.invm_min, kincuts_.invm_max );
cuts_.cuts.central.pt_sum.save( (bool&)kincuts_.iptsum, kincuts_.ptsum_min, kincuts_.ptsum_max );
cuts_.cuts.central.rapidity_diff.save( (bool&)kincuts_.idely, kincuts_.dely_min, kincuts_.dely_max );
//===========================================================================================
// feed run parameters to the common block
//===========================================================================================
params_.icontri = (int)cuts_.mode;
params_.imethod = method_;
params_.sfmod = (int)cuts_.structure_functions->type;
params_.pdg_l = pair_;
//-------------------------------------------------------------------------------------------
// incoming beams information
//-------------------------------------------------------------------------------------------
params_.inp1 = cuts_.incoming_beams.first.pz;
params_.inp2 = cuts_.incoming_beams.second.pz;
const HeavyIon in1 = (HeavyIon)cuts_.incoming_beams.first.pdg;
if ( in1 ) {
params_.a_nuc1 = in1.A;
params_.z_nuc1 = (unsigned short)in1.Z;
if ( params_.z_nuc1 > 1 ) {
event_->getOneByRole( Particle::IncomingBeam1 ).setPdgId( (PDG)in1 );
event_->getOneByRole( Particle::OutgoingBeam1 ).setPdgId( (PDG)in1 );
}
}
else
params_.a_nuc1 = params_.z_nuc1 = 1;
const HeavyIon in2 = (HeavyIon)cuts_.incoming_beams.second.pdg;
if ( in2 ) {
params_.a_nuc2 = in2.A;
params_.z_nuc2 = (unsigned short)in2.Z;
if ( params_.z_nuc2 > 1 ) {
event_->getOneByRole( Particle::IncomingBeam2 ).setPdgId( (PDG)in2 );
event_->getOneByRole( Particle::OutgoingBeam2 ).setPdgId( (PDG)in2 );
}
}
else
params_.a_nuc2 = params_.z_nuc2 = 1;
//-------------------------------------------------------------------------------------------
// intermediate partons information
//-------------------------------------------------------------------------------------------
params_.iflux1 = (int)cuts_.incoming_beams.first.kt_flux;
params_.iflux2 = (int)cuts_.incoming_beams.second.kt_flux;
if ( (KTFlux)params_.iflux1 == KTFlux::P_Gluon_KMR )
event_->getOneByRole( Particle::Parton1 ).setPdgId( PDG::gluon );
if ( (KTFlux)params_.iflux2 == KTFlux::P_Gluon_KMR )
event_->getOneByRole( Particle::Parton2 ).setPdgId( PDG::gluon );
}
double
FortranKTProcess::computeKTFactorisedMatrixElement()
{
ktkin_.q1t = qt1_;
ktkin_.q2t = qt2_;
ktkin_.phiq1t = phi_qt1_;
ktkin_.phiq2t = phi_qt2_;
ktkin_.y1 = y1_;
ktkin_.y2 = y2_;
ktkin_.ptdiff = pt_diff_;
ktkin_.phiptdiff = phi_pt_diff_;
ktkin_.m_x = MX_;
ktkin_.m_y = MY_;
//--- compute the event weight
double weight = 0.;
func_( weight );
return weight;
}
void
FortranKTProcess::fillCentralParticlesKinematics()
{
//===========================================================================================
// outgoing beam remnants
//===========================================================================================
PX_ = Particle::Momentum( evtkin_.px );
PY_ = Particle::Momentum( evtkin_.py );
// express these momenta per nucleon
PX_ *= 1./params_.a_nuc1;
PY_ *= 1./params_.a_nuc2;
//===========================================================================================
// intermediate partons
//===========================================================================================
const Particle::Momentum mom_par1 = mom_ip1_-PX_, mom_par2 = mom_ip2_-PY_;
event_->getOneByRole( Particle::Parton1 ).setMomentum( mom_par1 );
event_->getOneByRole( Particle::Parton2 ).setMomentum( mom_par2 );
event_->getOneByRole( Particle::Intermediate ).setMomentum( mom_par1+mom_par2 );
//===========================================================================================
// central system
//===========================================================================================
Particles& oc = event_->getByRole( Particle::CentralSystem );
for ( int i = 0; i < evtkin_.nout; ++i ) {
Particle& p = oc[i];
p.setPdgId( evtkin_.pdg[i] );
p.setStatus( Particle::Status::FinalState );
p.setMomentum( Particle::Momentum( evtkin_.pc[i] ) );
}
}
}
}
-
diff --git a/CepGen/Processes/FortranKTProcess.h b/CepGen/Processes/FortranKTProcess.h
index b36674a..f00bd37 100644
--- a/CepGen/Processes/FortranKTProcess.h
+++ b/CepGen/Processes/FortranKTProcess.h
@@ -1,38 +1,37 @@
#ifndef CepGen_Processes_FortranKTProcess_h
#define CepGen_Processes_FortranKTProcess_h
#include "CepGen/Processes/GenericKTProcess.h"
#include <functional>
namespace CepGen
{
- namespace Process
+ namespace process
{
/// Compute the matrix element for a generic \f$k_T\f$-factorised process defined in a Fortran subroutine
class FortranKTProcess : public GenericKTProcess
{
public:
FortranKTProcess( const ParametersList& params, const char* name, const char* descr, std::function<void(double&)> func );
ProcessPtr clone( const ParametersList& params ) const override { return ProcessPtr( new FortranKTProcess( *this ) ); }
private:
void preparePhaseSpace() override;
double computeKTFactorisedMatrixElement() override;
void fillCentralParticlesKinematics() override;
int pair_; ///< Outgoing particles type
int method_; ///< Computation method for the process
std::function<void(double&)> func_; ///< Subroutine to be called for weight computation
double y1_; ///< First outgoing particle rapidity
double y2_; ///< Second outgoing particle rapidity
double pt_diff_; ///< Transverse momentum balance between outgoing particles
double phi_pt_diff_; ///< Azimutal angle difference between outgoing particles
Particle::Momentum mom_ip1_; ///< First incoming beam momentum
Particle::Momentum mom_ip2_; ///< Second incoming beam momentum
};
}
}
#endif
-
diff --git a/CepGen/Processes/FortranProcesses.h b/CepGen/Processes/FortranProcesses.h
index ec832b8..54e0db5 100644
--- a/CepGen/Processes/FortranProcesses.h
+++ b/CepGen/Processes/FortranProcesses.h
@@ -1,18 +1,17 @@
#ifndef CepGen_Processes_FortranProcesses_h
#define CepGen_Processes_FortranProcesses_h
#include "CepGen/Processes/FortranKTProcess.h"
#include "CepGen/Processes/ProcessesHandler.h"
#include "CepGen/Core/ParametersList.h"
#define DECLARE_FORTRAN_SUBROUTINE( method ) \
extern "C" { extern void method ## _( double& ); }
#define PROCESS_F77_NAME( name ) F77_ ## name
#define REGISTER_FORTRAN_PROCESS( name, method, description ) \
- struct PROCESS_F77_NAME( name ) : public CepGen::Process::FortranKTProcess { \
- PROCESS_F77_NAME( name )() : CepGen::Process::FortranKTProcess( CepGen::ParametersList(), STRINGIFY( name ), description, method ## _ ) {} }; \
+ struct PROCESS_F77_NAME( name ) : public CepGen::process::FortranKTProcess { \
+ PROCESS_F77_NAME( name )() : CepGen::process::FortranKTProcess( CepGen::ParametersList(), STRINGIFY( name ), description, method ## _ ) {} }; \
REGISTER_PROCESS( name, PROCESS_F77_NAME( name ) )
#endif
-
diff --git a/CepGen/Processes/GamGamLL.cpp b/CepGen/Processes/GamGamLL.cpp
index 34ad076..52c0661 100644
--- a/CepGen/Processes/GamGamLL.cpp
+++ b/CepGen/Processes/GamGamLL.cpp
@@ -1,1150 +1,1149 @@
#include "CepGen/Processes/GamGamLL.h"
#include "CepGen/Event/Event.h"
#include "CepGen/Core/Exception.h"
#include "CepGen/Physics/Constants.h"
#include "CepGen/Physics/FormFactors.h"
#include "CepGen/Physics/PDG.h"
#include "CepGen/Processes/ProcessesHandler.h"
namespace CepGen
{
- namespace Process
+ namespace process
{
//---------------------------------------------------------------------------------------------
GamGamLL::Masses::Masses() :
MX2_( 0. ), MY2_( 0. ), Ml2_( 0. ),
w12_( 0. ), w31_( 0. ), dw31_( 0. ), w52_( 0. ), dw52_( 0. )
{}
//---------------------------------------------------------------------------------------------
GamGamLL::GamGamLL( const ParametersList& params ) :
GenericProcess( "lpair", "pp → p(*) ( ɣɣ → l⁺l¯ ) p(*)" ),
n_opt_( params.get<int>( "nopt", 0 ) ),
pair_( params.get<int>( "pair", 0 ) ),
ep1_( 0. ), ep2_( 0. ), p_cm_( 0. ),
ec4_( 0. ), pc4_( 0. ), mc4_( 0. ), w4_( 0. ),
p12_( 0. ), p1k2_( 0. ), p2k1_( 0. ),
p13_( 0. ), p14_( 0. ), p25_( 0. ),
q1dq_( 0. ), q1dq2_( 0. ),
s1_( 0. ), s2_( 0. ),
epsi_( 0. ),
g5_( 0. ), g6_( 0. ), a5_( 0. ), a6_( 0. ), bb_( 0. ),
gram_( 0. ),
dd1_( 0. ), dd2_( 0. ), dd3_( 0. ), dd4_( 0. ), dd5_( 0. ),
delta_( 0. ),
g4_( 0. ), sa1_( 0. ), sa2_( 0. ),
sl1_( 0. ),
cos_theta4_( 0. ), sin_theta4_( 0. ),
al4_( 0. ), be4_( 0. ), de3_( 0. ), de5_( 0. ),
pt4_( 0. ),
jacobian_( 0. )
{}
//---------------------------------------------------------------------------------------------
void
GamGamLL::addEventContent()
{
GenericProcess::setEventContent( {
{ Particle::IncomingBeam1, PDG::proton },
{ Particle::IncomingBeam2, PDG::proton },
{ Particle::Parton1, PDG::photon },
{ Particle::Parton2, PDG::photon }
}, {
{ Particle::OutgoingBeam1, { PDG::proton } },
{ Particle::OutgoingBeam2, { PDG::proton } },
{ Particle::CentralSystem, { (PDG)pair_, (PDG)pair_ } }
} );
}
unsigned int
GamGamLL::numDimensions() const
{
switch ( cuts_.mode ) {
case KinematicsMode::ElectronProton: default:
throw CG_FATAL( "GamGamLL" )
<< "Process mode " << cuts_.mode << " not (yet) supported! "
<< "Please contact the developers to consider an implementation.";
case KinematicsMode::ElasticElastic:
return 7;
case KinematicsMode::ElasticInelastic:
case KinematicsMode::InelasticElastic:
return 8;
case KinematicsMode::InelasticInelastic:
return 9;
}
}
//---------------------------------------------------------------------------------------------
void
GamGamLL::setKinematics( const Kinematics& kin )
{
GenericProcess::setKinematics( kin );
masses_.Ml2_ = event_->getByRole( Particle::CentralSystem )[0].mass2();
w_limits_ = cuts_.cuts.central.mass_single;
if ( !w_limits_.hasMax() )
w_limits_.max() = s_;
// The minimal energy for the central system is its outgoing leptons' mass energy (or wmin_ if specified)
if ( !w_limits_.hasMin() )
w_limits_.min() = 4.*masses_.Ml2_;
// The maximal energy for the central system is its CM energy with the outgoing particles' mass energy substracted (or wmax if specified)
w_limits_.max() = std::min( pow( sqs_-MX_-MY_, 2 ), w_limits_.max() );
CG_DEBUG_LOOP( "GamGamLL:setKinematics" )
<< "w limits = " << w_limits_ << "\n\t"
<< "wmax/wmin = " << w_limits_.max()/w_limits_.min();
q2_limits_ = cuts_.cuts.initial.q2;
mx_limits_ = cuts_.cuts.remnants.mass_single;
}
//---------------------------------------------------------------------------------------------
bool
GamGamLL::pickin()
{
CG_DEBUG_LOOP( "GamGamLL" )
<< "Optimised mode? " << n_opt_;
jacobian_ = 0.;
w4_ = mc4_*mc4_;
// sig1 = sigma and sig2 = sigma' in [1]
const double sig = mc4_+MY_;
double sig1 = sig*sig,
sig2 = sig1;
CG_DEBUG_LOOP( "GamGamLL" )
<< "mc4 = " << mc4_ << "\n\t"
<< "sig1 = " << sig1 << "\n\t"
<< "sig2 = " << sig2;
// Mass difference between the first outgoing particle
// and the first incoming particle
masses_.w31_ = masses_.MX2_-w1_;
// Mass difference between the second outgoing particle
// and the second incoming particle
masses_.w52_ = masses_.MY2_-w2_;
// Mass difference between the two incoming particles
masses_.w12_ = w1_-w2_;
// Mass difference between the central two-photons system
// and the second outgoing particle
const double d6 = w4_-masses_.MY2_;
CG_DEBUG_LOOP( "GamGamLL" )
<< "w1 = " << w1_ << "\n\t"
<< "w2 = " << w2_ << "\n\t"
<< "w3 = " << masses_.MX2_ << "\n\t"
<< "w4 = " << w4_ << "\n\t"
<< "w5 = " << masses_.MY2_;;
CG_DEBUG_LOOP( "GamGamLL" )
<< "w31 = " << masses_.w31_ << "\n\t"
<< "w52 = " << masses_.w52_ << "\n\t"
<< "w12 = " << masses_.w12_;;
const double ss = s_+masses_.w12_;
const double rl1 = ss*ss-4.*w1_*s_; // lambda(s, m1**2, m2**2)
if ( rl1 <= 0. ) {
CG_WARNING( "GamGamLL" ) << "rl1 = " << rl1 << " <= 0";
return false;
}
sl1_ = sqrt( rl1 );
s2_ = 0.;
double ds2 = 0.;
if ( n_opt_ == 0 ) {
const double smax = s_+masses_.MX2_-2.*MX_*sqs_;
map( x(2), Limits( sig1, smax ), s2_, ds2, "s2" );
sig1 = s2_; //FIXME!!!!!!!!!!!!!!!!!!!!
}
CG_DEBUG_LOOP( "GamGamLL" )
<< "s2 = " << s2_;
const double sp = s_+masses_.MX2_-sig1,
d3 = sig1-w2_;
const double rl2 = sp*sp-4.*s_*masses_.MX2_; // lambda(s, m3**2, sigma)
if ( rl2 <= 0. ) {
CG_DEBUG( "GamGamLL" ) << "rl2 = " << rl2 << " <= 0";
return false;
}
const double sl2 = sqrt( rl2 );
double t1_max = w1_+masses_.MX2_-( ss*sp+sl1_*sl2 )/( 2.*s_ ), // definition from eq. (A.4) in [1]
t1_min = ( masses_.w31_*d3+( d3-masses_.w31_ )*( d3*w1_-masses_.w31_*w2_ )/s_ )/t1_max; // definition from eq. (A.5) in [1]
// FIXME dropped in CDF version
if ( t1_max > -q2_limits_.min() ) {
CG_WARNING( "GamGamLL" ) << "t1max = " << t1_max << " > -q2min = " << ( -q2_limits_.min() );
return false;
}
if ( t1_min < -q2_limits_.max() && q2_limits_.hasMax() ) {
CG_DEBUG( "GamGamLL" ) << "t1min = " << t1_min << " < -q2max = " << -q2_limits_.max();
return false;
}
if ( t1_max < -q2_limits_.max() && q2_limits_.hasMax() )
t1_max = -q2_limits_.max();
if ( t1_min > -q2_limits_.min() && q2_limits_.hasMin() )
t1_min = -q2_limits_.min();
/////
// t1, the first photon propagator, is defined here
t1_ = 0.;
double dt1 = 0.;
map( x(0), Limits( t1_min, t1_max ), t1_, dt1, "t1" );
// changes wrt mapt1 : dx->-dx
dt1 *= -1.;
CG_DEBUG_LOOP( "GamGamLL" )
<< "Definition of t1 = " << t1_ << " according to\n\t"
<< "(t1min, t1max) = (" << t1_min << ", " << t1_max << ")";
dd4_ = w4_-t1_;
const double d8 = t1_-w2_, t13 = t1_-w1_-masses_.MX2_;
sa1_ = -pow( t1_-masses_.w31_, 2 )/4.+w1_*t1_;
if ( sa1_ >= 0. ) {
CG_WARNING( "GamGamLL" ) << "sa1_ = " << sa1_ << " >= 0";
return false;
}
const double sl3 = sqrt( -sa1_ );
// one computes splus and (s2x=s2max)
double splus, s2max;
if ( w1_ != 0. ) {
const double inv_w1 = 1./w1_;
const double sb = masses_.MX2_ + 0.5 * ( s_*( t1_-masses_.w31_ )+masses_.w12_*t13 )*inv_w1,
sd = sl1_*sl3*inv_w1,
se =( s_*( t1_*( s_+t13-w2_ )-w2_*masses_.w31_ )+masses_.MX2_*( masses_.w12_*d8+w2_*masses_.MX2_ ) )*inv_w1;
if ( fabs( ( sb-sd )/sd ) >= 1. ) {
splus = sb-sd;
s2max = se/splus;
}
else {
s2max = sb+sd;
splus = se/s2max;
}
}
else { // 3
s2max = ( s_*( t1_*( s_+d8-masses_.MX2_ )-w2_*masses_.MX2_ )+w2_*masses_.MX2_*( w2_+masses_.MX2_-t1_ ) )/( ss*t13 );
splus = sig2;
}
// 4
double s2x = s2max;
CG_DEBUG_LOOP( "GamGamLL" )
<< "s2x = s2max = " << s2x;
if ( n_opt_ < 0 ) { // 5
if ( splus > sig2 ) {
sig2 = splus;
CG_DEBUG_LOOP( "GamGamLL" )
<< "sig2 truncated to splus = " << splus;
}
if ( n_opt_ < -1 )
map( x(2), Limits( sig2, s2max ), s2_, ds2, "s2" );
else
mapla( t1_, w2_, x(2), sig2, s2max, s2_, ds2 ); // n_opt_==-1
s2x = s2_;
}
else if ( n_opt_ == 0 )
s2x = s2_; // 6
CG_DEBUG_LOOP( "GamGamLL" )
<< "s2x = " << s2x;
// 7
const double r1 = s2x-d8, r2 = s2x-d6;
const double rl4 = ( r1*r1-4.*w2_*s2x )*( r2*r2-4.*masses_.MY2_*s2x );
if ( rl4 <= 0. ) {
CG_DEBUG_LOOP( "GamGamLL" )
<< "rl4 = " << rl4 << " <= 0";
return false;
}
const double sl4 = sqrt( rl4 );
// t2max, t2min definitions from eq. (A.12) and (A.13) in [1]
const double t2_max = w2_+masses_.MY2_-( r1*r2+sl4 )/s2x * 0.5,
t2_min = ( masses_.w52_*dd4_+( dd4_-masses_.w52_ )*( dd4_*w2_-masses_.w52_*t1_ )/s2x )/t2_max;
// t2, the second photon propagator, is defined here
t2_ = 0.;
double dt2 = 0.;
map( x(1), Limits( t2_min, t2_max ), t2_, dt2, "t2" );
// changes wrt mapt2 : dx->-dx
dt2 *= -1.;
// \f$\delta_6=m_4^2-m_5^2\f$ as defined in Vermaseren's paper
const double tau = t1_-t2_,
r3 = dd4_-t2_,
r4 = masses_.w52_-t2_;
CG_DEBUG_LOOP( "GamGamLL" )
<< "r1 = " << r1 << "\n\t"
<< "r2 = " << r2 << "\n\t"
<< "r3 = " << r3 << "\n\t"
<< "r4 = " << r4;
const double b = r3*r4-2.*( t1_+w2_ )*t2_;
const double c = t2_*d6*d8+( d6-d8 )*( d6*w2_-d8*masses_.MY2_ );
const double t25 = t2_-w2_-masses_.MY2_;
sa2_ = -0.25 * r4*r4 + w2_*t2_;
if ( sa2_ >= 0. ) {
CG_WARNING( "GamGamLL" )
<< "sa2_ = " << sa2_ << " >= 0";
return false;
}
const double sl6 = 2.*sqrt( -sa2_ );
g4_ = -r3*r3/4.+t1_*t2_;
if ( g4_ >= 0. ) {
CG_WARNING( "GamGamLL" ) << "g4_ = " << g4_ << " >= 0";
return false;
}
const double sl7 = 2.*sqrt( -g4_ ),
sl5 = sl6*sl7;
double s2p, s2min;
if ( fabs( ( sl5-b )/sl5 ) >= 1. ) {
s2p = 0.5 * ( sl5-b )/t2_;
s2min = c/( t2_*s2p );
}
else { // 8
s2min = 0.5 * ( -sl5-b )/t2_;
s2p = c/( t2_*s2min );
}
// 9
if ( n_opt_ > 1 )
map( x( 2 ), Limits( s2min, s2max ), s2_, ds2, "s2" );
else if ( n_opt_ == 1 )
mapla( t1_, w2_, x( 2 ), s2min, s2max, s2_, ds2 );
const double ap = -0.25*pow( s2_+d8, 2 )+s2_*t1_;
CG_DEBUG_LOOP( "GamGamLL" )
<< "s2 = " << s2_ << ", s2max = " << s2max << ", splus = " << splus;
if ( w1_ != 0. )
dd1_ = -0.25 * ( s2_-s2max ) * ( s2_-splus ) * w1_; // 10
else
dd1_ = 0.25 * ( s2_-s2max ) * ss * t13;
// 11
dd2_ = -0.25 * t2_*( s2_-s2p )*( s2_-s2min );
CG_DEBUG_LOOP( "GamGamLL" )
<< "t2 = " << t2_ << "\n\t"
<< "s2 = " << s2_ << "\n\t"
<< "s2p = " << s2p << "\n\t"
<< "s2min = " << s2min << "\n\t"
<< "dd2 = " << dd2_;;
const double yy4 = cos( M_PI*x( 3 ) );
const double dd = dd1_*dd2_;
p12_ = 0.5 * ( s_-w1_-w2_ );
const double st = s2_-t1_-w2_;
const double delb = ( 2.*w2_*r3+r4*st )*( 4.*p12_*t1_-( t1_-masses_.w31_ )*st )/( 16.*ap );
if ( dd <= 0. ) {
CG_DEBUG_LOOP( "GamGamLL" )
<< std::scientific
<< "dd = " << dd << " <= 0\n\t"
<< "dd1 = " << dd1_ << "\t"
<< "dd2 = " << dd2_ << std::fixed;
return false;
}
delta_ = delb - yy4*st*sqrt( dd )/ ap * 0.5;
s1_ = t2_+w1_+( 2.*p12_*r3-4.*delta_ )/st;
if ( ap >= 0. ) {
CG_DEBUG_LOOP( "GamGamLL" )
<< "ap = " << ap << " >= 0";
return false;
}
jacobian_ = ds2 * dt1 * dt2 * 0.125 * M_PI*M_PI/( sl1_*sqrt( -ap ) );
CG_DEBUG_LOOP( "GamGamLL" )
<< "Jacobian = " << std::scientific << jacobian_ << std::fixed;
gram_ = ( 1.-yy4*yy4 )*dd/ap;
p13_ = -0.5 * t13;
p14_ = 0.5 * ( tau+s1_-masses_.MX2_ );
p25_ = -0.5 * t25;
p1k2_ = 0.5 * ( s1_-t2_-w1_ );
p2k1_ = 0.5 * st;
if ( w2_ != 0. ) {
const double inv_w2 = 1./w2_;
const double sbb = 0.5 * ( s_*( t2_-masses_.w52_ )-masses_.w12_*t25 )*inv_w2 + masses_.MY2_,
sdd = 0.5 * sl1_*sl6*inv_w2,
see = ( s_*( t2_*( s_+t25-w1_ )-w1_*masses_.w52_ )+masses_.MY2_*( w1_*masses_.MY2_-masses_.w12_*( t2_-w1_ ) ) )*inv_w2;
double s1m = 0., s1p = 0.;
if ( sbb/sdd >= 0. ) {
s1p = sbb+sdd;
s1m = see/s1p;
}
else {
s1m = sbb-sdd;
s1p = see/s1m;
} // 12
dd3_ = -0.25 * w2_*( s1p-s1_ )*( s1m-s1_ ); // 13
}
else { // 14
const double s1p = ( s_*( t2_*( s_-masses_.MY2_+t2_-w1_ )-w1_*masses_.MY2_ )+w1_*masses_.MY2_*( w1_+masses_.MY2_-t2_ ) )/( t25*( s_-masses_.w12_ ) );
dd3_ = -0.25 * t25*( s_-masses_.w12_ )*( s1p-s1_ );
}
// 15
//const double acc3 = (s1p-s1_)/(s1p+s1_);
const double ssb = t2_+0.5 * w1_-r3*( masses_.w31_-t1_ )/t1_,
ssd = sl3*sl7/t1_,
sse = ( t2_-w1_ )*( w4_-masses_.MX2_ )+( t2_-w4_+masses_.w31_ )*( ( t2_-w1_)*masses_.MX2_-(w4_-masses_.MX2_)*w1_)/t1_;
double s1pp, s1pm;
if ( ssb/ssd >= 0. ) {
s1pp = ssb+ssd;
s1pm = sse/s1pp;
}
else { // 16
s1pm = ssb-ssd;
s1pp = sse/s1pm;
}
// 17
dd4_ = -0.25 * t1_*( s1_-s1pp )*( s1_-s1pm );
//const double acc4 = ( s1_-s1pm )/( s1_+s1pm );
dd5_ = dd1_+dd3_+( ( p12_*( t1_-masses_.w31_ )*0.5-w1_*p2k1_ )*( p2k1_*( t2_-masses_.w52_ )-w2_*r3 )
-delta_*( 2.*p12_*p2k1_-w2_*( t1_-masses_.w31_ ) ) ) / p2k1_;
return true;
}
//---------------------------------------------------------------------------------------------
bool
GamGamLL::orient()
{
if ( !pickin() || jacobian_ == 0. ) {
CG_DEBUG_LOOP( "GamGamLL" )
<< "Pickin failed! Jacobian = " << jacobian_;
return false;
}
const double re = 0.5 / sqs_;
ep1_ = re*( s_+masses_.w12_ );
ep2_ = re*( s_-masses_.w12_ );
CG_DEBUG_LOOP( "GamGamLL" )
<< std::scientific
<< " re = " << re << "\n\t"
<< "w12_ = " << masses_.w12_
<< std::fixed;
CG_DEBUG_LOOP( "GamGamLL" )
<< "Incoming particles' energy = " << ep1_ << ", " << ep2_;
p_cm_ = re*sl1_;
de3_ = re*(s2_-masses_.MX2_+masses_.w12_);
de5_ = re*(s1_-masses_.MY2_-masses_.w12_);
// Final state energies
const double ep3 = ep1_-de3_,
ep5 = ep2_-de5_;
ec4_ = de3_+de5_;
if ( ec4_ < mc4_ ) {
CG_WARNING( "GamGamLL" )
<< "ec4_ = " << ec4_ << " < mc4_ = " << mc4_ << "\n\t"
<< "==> de3 = " << de3_ << ", de5 = " << de5_;
return false;
}
// What if the protons' momenta are not along the z-axis?
pc4_ = sqrt( ec4_*ec4_-mc4_*mc4_ );
if ( pc4_ == 0. ) {
CG_WARNING( "GamGamLL" ) << "pzc4 is null and should not be...";
return false;
}
const double pp3 = sqrt( ep3*ep3-masses_.MX2_ ), pt3 = sqrt( dd1_/s_ )/p_cm_,
pp5 = sqrt( ep5*ep5-masses_.MY2_ ), pt5 = sqrt( dd3_/s_ )/p_cm_;
CG_DEBUG_LOOP( "GamGamLL" )
<< "Central system's energy: E4 = " << ec4_ << "\n\t"
<< " momentum: p4 = " << pc4_ << "\n\t"
<< " invariant mass: m4 = " << mc4_ << "\n\t"
<< "Outgoing particles' energy: E3 = " << ep3 << "\n\t"
<< " E5 = " << ep5;
const double sin_theta3 = pt3/pp3,
sin_theta5 = pt5/pp5;
CG_DEBUG_LOOP( "GamGamLL" )
<< std::scientific
<< "sin_theta3 = " << sin_theta3 << "\n\t"
<< "sin_theta5 = " << sin_theta5
<< std::fixed;
if ( sin_theta3 > 1. ) {
CG_WARNING( "GamGamLL" )
<< "sin_theta3 = " << sin_theta3 << " > 1";
return false;
}
if ( sin_theta5 > 1. ) {
CG_WARNING( "GamGamLL" )
<< "sin_theta5 = " << sin_theta5 << " > 1";
return false;
}
double ct3 = sqrt( 1.-sin_theta3*sin_theta3 ),
ct5 = sqrt( 1.-sin_theta5*sin_theta5 );
if ( ep1_*ep3 < p13_ )
ct3 *= -1.;
if ( ep2_*ep5 > p25_ )
ct5 *= -1.;
CG_DEBUG_LOOP( "GamGamLL" )
<< "ct3 = " << ct3 << "\n\t"
<< "ct5 = " << ct5;
if ( dd5_ < 0. ) {
CG_WARNING( "GamGamLL" )
<< "dd5 = " << dd5_ << " < 0";
return false;
}
// Centre of mass system kinematics (theta4 and phi4)
pt4_ = sqrt( dd5_/s_ )/p_cm_;
sin_theta4_ = pt4_/pc4_;
if ( sin_theta4_ > 1. ) {
CG_WARNING( "GamGamLL" )
<< "st4 = " << sin_theta4_ << " > 1";
return false;
}
cos_theta4_ = sqrt( 1.-sin_theta4_*sin_theta4_ );
if ( ep1_*ec4_ < p14_ )
cos_theta4_ *= -1.;
al4_ = 1.-cos_theta4_;
be4_ = 1.+cos_theta4_;
if ( cos_theta4_ < 0. ) be4_ = sin_theta4_*sin_theta4_/al4_;
else al4_ = sin_theta4_*sin_theta4_/be4_;
CG_DEBUG_LOOP( "GamGamLL" )
<< "ct4 = " << cos_theta4_ << "\n\t"
<< "al4 = " << al4_ << ", be4 = " << be4_;
const double rr = sqrt( -gram_/s_ )/( p_cm_*pt4_ );
const double sin_phi3 = rr / pt3,
sin_phi5 = -rr / pt5;
if ( fabs( sin_phi3 ) > 1. ) {
CG_WARNING( "GamGamLL" )
<< "sin(phi_3) = " << sin_phi3 << " while it must be in [-1 ; 1]";
return false;
}
if ( fabs( sin_phi5 ) > 1. ) {
CG_WARNING( "GamGamLL" )
<< "sin(phi_5) = " << sin_phi5 << " while it must be in [-1 ; 1]";
return false;
}
const double cos_phi3 = -sqrt( 1.-sin_phi3*sin_phi3 ),
cos_phi5 = -sqrt( 1.-sin_phi5*sin_phi5 );
p3_lab_ = Particle::Momentum( pp3*sin_theta3*cos_phi3, pp3*sin_theta3*sin_phi3, pp3*ct3, ep3 );
p5_lab_ = Particle::Momentum( pp5*sin_theta5*cos_phi5, pp5*sin_theta5*sin_phi5, pp5*ct5, ep5 );
const double a1 = p3_lab_.px()-p5_lab_.px();
CG_DEBUG_LOOP( "GamGamLL" )
<< "Kinematic quantities\n\t"
<< "cos(theta3) = " << ct3 << "\t" << "sin(theta3) = " << sin_theta3 << "\n\t"
<< "cos( phi3 ) = " << cos_phi3 << "\t" << "sin( phi3 ) = " << sin_phi3 << "\n\t"
<< "cos(theta4) = " << cos_theta4_ << "\t" << "sin(theta4) = " << sin_theta4_ << "\n\t"
<< "cos(theta5) = " << ct5 << "\t" << "sin(theta5) = " << sin_theta5 << "\n\t"
<< "cos( phi5 ) = " << cos_phi5 << "\t" << "sin( phi5 ) = " << sin_phi5 << "\n\t"
<< "a1 = " << a1;
if ( fabs( pt4_+p3_lab_.px()+p5_lab_.px() ) < fabs( fabs( a1 )-pt4_ ) ) {
CG_DEBUG_LOOP( "GamGamLL" )
<< "|pt4+pt3*cos(phi3)+pt5*cos(phi5)| < | |a1|-pt4 |\n\t"
<< "pt4 = " << pt4_ << "\t"
<< "pt5 = " << pt5 << "\n\t"
<< "cos(phi3) = " << cos_phi3 << "\t"
<< "cos(phi5) = " << cos_phi5 << "\n\t"
<< "a1 = " << a1;
return true;
}
if ( a1 < 0. ) p5_lab_[0] = -p5_lab_.px();
else p3_lab_[0] = -p3_lab_.px();
return true;
}
//---------------------------------------------------------------------------------------------
double
GamGamLL::computeOutgoingPrimaryParticlesMasses( double x, double outmass, double lepmass, double& dw )
{
const double mx0 = mp_+ParticleProperties::mass( PDG::piZero ); // 1.07
const double wx2min = pow( std::max( mx0, mx_limits_.min() ), 2 ),
wx2max = pow( std::min( sqs_-outmass-2.*lepmass, mx_limits_.max() ), 2 );
double mx2 = 0., dmx2 = 0.;
map( x, Limits( wx2min, wx2max ), mx2, dmx2, "mx2" );
CG_DEBUG_LOOP( "GamGamLL" )
<< "mX^2 in range (" << wx2min << ", " << wx2max << "), x = " << x << "\n\t"
<< "mX^2 = " << mx2 << ", d(mX^2) = " << dmx2 << "\n\t"
<< "mX = " << sqrt( mx2 ) << ", d(mX) = " << sqrt( dmx2 );
dw = sqrt( dmx2 );
return sqrt( mx2 );
}
//---------------------------------------------------------------------------------------------
void
GamGamLL::beforeComputeWeight()
{
if ( !GenericProcess::is_point_set_ ) return;
const Particle& p1 = event_->getOneByRole( Particle::IncomingBeam1 ),
&p2 = event_->getOneByRole( Particle::IncomingBeam2 );
ep1_ = p1.energy();
ep2_ = p2.energy();
//std::cout << __PRETTY_FUNCTION__ << ":" << w_limits_ << "|" << q2_limits_ << "|" << mx_limits_ << std::endl;
switch ( cuts_.mode ) {
case KinematicsMode::ElectronProton: default:
CG_ERROR( "GamGamLL" ) << "Case not yet supported!"; break;
case KinematicsMode::ElasticElastic:
masses_.dw31_ = masses_.dw52_ = 0.; break;
case KinematicsMode::InelasticElastic: {
const double m = computeOutgoingPrimaryParticlesMasses( x( 7 ), p1.mass(), sqrt( masses_.Ml2_ ), masses_.dw31_ );
event_->getOneByRole( Particle::OutgoingBeam1 ).setMass( m );
event_->getOneByRole( Particle::OutgoingBeam2 ).setMass( ParticleProperties::mass( p2.pdgId() ) );
} break;
case KinematicsMode::ElasticInelastic: {
const double m = computeOutgoingPrimaryParticlesMasses( x( 7 ), p2.mass(), sqrt( masses_.Ml2_ ), masses_.dw52_ );
event_->getOneByRole( Particle::OutgoingBeam1 ).setMass( ParticleProperties::mass( p1.pdgId() ) );
event_->getOneByRole( Particle::OutgoingBeam2 ).setMass( m );
} break;
case KinematicsMode::InelasticInelastic: {
const double mx = computeOutgoingPrimaryParticlesMasses( x( 7 ), p2.mass(), sqrt( masses_.Ml2_ ), masses_.dw31_ );
event_->getOneByRole( Particle::OutgoingBeam1 ).setMass( mx );
const double my = computeOutgoingPrimaryParticlesMasses( x( 8 ), p1.mass(), sqrt( masses_.Ml2_ ), masses_.dw52_ );
event_->getOneByRole( Particle::OutgoingBeam2 ).setMass( my );
} break;
}
MX_ = event_->getOneByRole( Particle::OutgoingBeam1 ).mass();
MY_ = event_->getOneByRole( Particle::OutgoingBeam2 ).mass();
masses_.MX2_ = MX_*MX_;
masses_.MY2_ = MY_*MY_;
}
//---------------------------------------------------------------------------------------------
double
GamGamLL::computeWeight()
{
CG_DEBUG_LOOP( "GamGamLL" )
<< "sqrt(s) = " << sqs_ << " GeV\n\t"
<< "m(X1) = " << MX_ << " GeV\t"
<< "m(X2) = " << MY_ << " GeV";
// compute the two-photon energy for this point
w4_ = 0.;
double dw4 = 0.;
map( x( 4 ), w_limits_, w4_, dw4, "w4" );
mc4_ = sqrt( w4_ );
CG_DEBUG_LOOP( "GamGamLL" )
<< "Computed value for w4 = " << w4_ << " → mc4 = " << mc4_;
if ( !orient() )
return 0.;
if ( jacobian_ == 0. ) {
CG_WARNING( "GamGamLL" ) << "dj = " << jacobian_;
return 0.;
}
if ( t1_ > 0. ) {
CG_WARNING( "GamGamLL" ) << "t1 = " << t1_ << " > 0";
return 0.;
}
if ( t2_ > 0. ) {
CG_WARNING( "GamGamLL" ) << "t2 = " << t2_ << " > 0";
return 0.;
}
const double ecm6 = w4_ / ( 2.*mc4_ ),
pp6cm = sqrt( ecm6*ecm6-masses_.Ml2_ );
jacobian_ *= dw4*pp6cm/( mc4_*Constants::sconstb*s_ );
// Let the most obscure part of this code begin...
const double e1mp1 = w1_ / ( ep1_+p_cm_ ),
e3mp3 = masses_.MX2_ / ( p3_lab_.energy()+p3_lab_.p() );
const double al3 = pow( sin( p3_lab_.theta() ), 2 )/( 1.+( p3_lab_.theta() ) );
// 2-photon system kinematics ?!
const double eg = ( w4_+t1_-t2_ )/( 2.*mc4_ );
double pg = sqrt( eg*eg-t1_ );
const double pgx = -p3_lab_.px()*cos_theta4_-sin_theta4_*( de3_-e1mp1 + e3mp3 + p3_lab_.p()*al3 ),
pgy = -p3_lab_.py(),
pgz = mc4_*de3_/( ec4_+pc4_ )-ec4_*de3_*al4_/mc4_-p3_lab_.px()*ec4_*sin_theta4_/mc4_+ec4_*cos_theta4_/mc4_*( p3_lab_.p()*al3+e3mp3-e1mp1 );
CG_DEBUG_LOOP( "GamGamLL" ) << "pg = " << Particle::Momentum( pgx, pgy, pgz );
const double pgp = sqrt( pgx*pgx + pgy*pgy ), // outgoing proton (3)'s transverse momentum
pgg = sqrt( pgp*pgp + pgz*pgz ); // outgoing proton (3)'s momentum
if ( pgg > pgp*0.9 && pgg > pg )
pg = pgg; //FIXME ???
// angles for the 2-photon system ?!
const double cpg = pgx/pgp, spg = pgy/pgp;
const double stg = pgp/pg;
const int theta_sign = ( pgz>0. ) ? 1 : -1;
const double ctg = theta_sign*sqrt( 1.-stg*stg );
double xx6 = x( 5 );
const double amap = 0.5 * ( w4_-t1_-t2_ ),
bmap = 0.5 * sqrt( ( pow( w4_-t1_-t2_, 2 )-4.*t1_*t2_ )*( 1.-4.*masses_.Ml2_/w4_ ) ),
ymap = ( amap+bmap )/( amap-bmap ),
beta = pow( ymap, 2.*xx6-1. );
xx6 = 0.5 * ( 1. + amap/bmap*( beta-1. )/( beta+1. ) );
xx6 = std::max( 0., std::min( xx6, 1. ) ); // xx6 in [0., 1.]
CG_DEBUG_LOOP( "GamGamLL" )
<< "amap = " << amap << "\n\t"
<< "bmap = " << bmap << "\n\t"
<< "ymap = " << ymap << "\n\t"
<< "beta = " << beta;
// 3D rotation of the first outgoing lepton wrt the CM system
const double theta6cm = acos( 1.-2.*xx6 );
// match the Jacobian
jacobian_ *= ( amap+bmap*cos( theta6cm ) );
jacobian_ *= ( amap-bmap*cos( theta6cm ) );
jacobian_ /= amap;
jacobian_ /= bmap;
jacobian_ *= log( ymap );
jacobian_ *= 0.5;
CG_DEBUG_LOOP( "GamGamLL" ) << "Jacobian = " << jacobian_;
CG_DEBUG_LOOP( "GamGamLL" )
<< "ctcm6 = " << cos( theta6cm ) << "\n\t"
<< "stcm6 = " << sin( theta6cm );
const double phi6cm = 2.*M_PI*x( 6 );
// First outgoing lepton's 3-momentum in the centre of mass system
Particle::Momentum p6cm = Particle::Momentum::fromPThetaPhi( pp6cm, theta6cm, phi6cm );
CG_DEBUG_LOOP( "GamGamLL" ) << "p3cm6 = " << p6cm;
const double h1 = stg*p6cm.pz()+ctg*p6cm.px();
const double pc6z = ctg*p6cm.pz()-stg*p6cm.px(), pc6x = cpg*h1-spg*p6cm.py();
const double qcx = 2.*pc6x, qcz = 2.*pc6z;
// qcy == QCY is never defined
const double el6 = ( ec4_*ecm6+pc4_*pc6z ) / mc4_,
h2 = ( ec4_*pc6z+pc4_*ecm6 ) / mc4_;
CG_DEBUG_LOOP( "GamGamLL" ) << "h1 = " << h1 << "\n\th2 = " << h2;
// first outgoing lepton's 3-momentum
const double p6x = cos_theta4_*pc6x+sin_theta4_*h2,
p6y = cpg*p6cm.py()+spg*h1,
p6z = cos_theta4_*h2-sin_theta4_*pc6x;
// first outgoing lepton's kinematics
p6_cm_ = Particle::Momentum( p6x, p6y, p6z, el6 );
CG_DEBUG_LOOP( "GamGamLL" ) << "p6(cm) = " << p6_cm_;
const double hq = ec4_*qcz/mc4_;
const Particle::Momentum qve(
cos_theta4_*qcx+sin_theta4_*hq,
2.*p6y,
cos_theta4_*hq-sin_theta4_*qcx,
pc4_*qcz/mc4_ // energy
);
// Available energy for the second lepton is the 2-photon system's energy with the first lepton's energy removed
const double el7 = ec4_-el6;
CG_DEBUG_LOOP( "GamGamLL" )
<< "Outgoing kinematics\n\t"
<< " first outgoing lepton: p = " << p6_cm_.p() << ", E = " << p6_cm_.energy() << "\n\t"
<< "second outgoing lepton: p = " << p7_cm_.p() << ", E = " << p7_cm_.energy();;
// Second outgoing lepton's 3-momentum
const double p7x = -p6x + pt4_,
p7y = -p6y,
p7z = -p6z + pc4_*cos_theta4_;
// second outgoing lepton's kinematics
p7_cm_ = Particle::Momentum( p7x, p7y, p7z, el7 );
//p6_cm_ = Particle::Momentum(pl6*st6*cp6, pl6*st6*sp6, pl6*ct6, el6);
//p7_cm_ = Particle::Momentum(pl7*st7*cp7, pl7*st7*sp7, pl7*ct7, el7);
q1dq_ = eg*( 2.*ecm6-mc4_ )-2.*pg*p6cm.pz();
q1dq2_ = ( w4_-t1_-t2_ ) * 0.5;
const double phi3 = p3_lab_.phi(), cos_phi3 = cos( phi3 ), sin_phi3 = sin( phi3 ),
phi5 = p5_lab_.phi(), cos_phi5 = cos( phi5 ), sin_phi5 = sin( phi5 );
bb_ = t1_*t2_+( w4_*pow( sin( theta6cm ), 2 ) + 4.*masses_.Ml2_*pow( cos( theta6cm ), 2 ) )*pg*pg;
const double c1 = p3_lab_.pt() * ( qve.px()*sin_phi3 - qve.py()*cos_phi3 ),
c2 = p3_lab_.pt() * ( qve.pz()*ep1_ - qve.energy() *p_cm_ ),
c3 = ( masses_.w31_*ep1_*ep1_ + 2.*w1_*de3_*ep1_ - w1_*de3_*de3_ + p3_lab_.pt2()*ep1_*ep1_ ) / ( p3_lab_.energy()*p_cm_ + p3_lab_.pz()*ep1_ );
const double b1 = p5_lab_.pt() * ( qve.px()*sin_phi5 - qve.py()*cos_phi5 ),
b2 = p5_lab_.pt() * ( qve.pz()*ep2_ + qve.energy() *p_cm_ ),
b3 = ( masses_.w52_*ep2_*ep2_ + 2.*w2_*de5_*ep2_ - w2_*de5_*de5_ + p5_lab_.pt2()*ep2_*ep2_ ) / ( ep2_*p5_lab_.pz() - p5_lab_.energy()*p_cm_ );
const double r12 = c2*sin_phi3 + qve.py()*c3,
r13 = -c2*cos_phi3 - qve.px()*c3;
const double r22 = b2*sin_phi5 + qve.py()*b3,
r23 = -b2*cos_phi5 - qve.px()*b3;
epsi_ = p12_*c1*b1 + r12*r22 + r13*r23;
g5_ = w1_*c1*c1 + r12*r12 + r13*r13;
g6_ = w2_*b1*b1 + r22*r22 + r23*r23;
const double pt3 = p3_lab_.pt(), pt5 = p5_lab_.pt();
a5_ = -( qve.px()*cos_phi3 + qve.py()*sin_phi3 )*pt3*p1k2_
-( ep1_*qve.energy()-p_cm_*qve.pz() )*( cos_phi3*cos_phi5 + sin_phi3*sin_phi5 )*pt3*pt5
+( de5_*qve.pz()+qve.energy()*( p_cm_+p5_lab_.pz() ) )*c3;
a6_ = -( qve.px()*cos_phi5 + qve.py()*sin_phi5 )*pt5*p2k1_
-( ep2_*qve.energy()+p_cm_*qve.pz() )*( cos_phi3*cos_phi5 + sin_phi3*sin_phi5 )*pt3*pt5
+( de3_*qve.pz()-qve.energy()*( p_cm_-p3_lab_.pz() ) )*b3;
CG_DEBUG_LOOP( "GamGamLL" )
<< "a5 = " << a5_ << "\n\t"
<< "a6 = " << a6_;
////////////////////////////////////////////////////////////////
// END of GAMGAMLL subroutine in the FORTRAN version
////////////////////////////////////////////////////////////////
const Particle::Momentum cm = event_->getOneByRole( Particle::IncomingBeam1 ).momentum()
+ event_->getOneByRole( Particle::IncomingBeam2 ).momentum();
////////////////////////////////////////////////////////////////
// INFO from f.f
////////////////////////////////////////////////////////////////
const double gamma = cm.energy() / sqs_, betgam = cm.pz() / sqs_;
//--- kinematics computation for both leptons
p6_cm_.betaGammaBoost( gamma, betgam );
p7_cm_.betaGammaBoost( gamma, betgam );
//--- cut on mass of final hadronic system (MX/Y)
if ( mx_limits_.valid() ) {
if ( ( cuts_.mode == KinematicsMode::InelasticElastic
|| cuts_.mode == KinematicsMode::InelasticInelastic )
&& !mx_limits_.passes( MX_ ) )
return 0.;
if ( ( cuts_.mode == KinematicsMode::ElasticInelastic
|| cuts_.mode == KinematicsMode::InelasticInelastic )
&& !mx_limits_.passes( MY_ ) )
return 0.;
}
//--- cut on the proton's Q2 (first photon propagator T1)
if ( !cuts_.cuts.initial.q2.passes( -t1_ ) )
return 0.;
//--- cuts on outgoing leptons' kinematics
if ( !cuts_.cuts.central.mass_sum.passes( ( p6_cm_+p7_cm_ ).mass() ) )
return 0.;
//----- cuts on the individual leptons
if ( cuts_.cuts.central.pt_single.valid() ) {
const Limits& pt_limits = cuts_.cuts.central.pt_single;
if ( !pt_limits.passes( p6_cm_.pt() ) || !pt_limits.passes( p7_cm_.pt() ) )
return 0.;
}
if ( cuts_.cuts.central.energy_single.valid() ) {
const Limits& energy_limits = cuts_.cuts.central.energy_single;
if ( !energy_limits.passes( p6_cm_.energy() ) || !energy_limits.passes( p7_cm_.energy() ) )
return 0.;
}
if ( cuts_.cuts.central.eta_single.valid() ) {
const Limits& eta_limits = cuts_.cuts.central.eta_single;
if ( !eta_limits.passes( p6_cm_.eta() ) || !eta_limits.passes( p7_cm_.eta() ) )
return 0.;
}
//--- compute the structure functions factors
switch ( cuts_.mode ) { // inherited from CDF version
case KinematicsMode::ElectronProton: default: jacobian_ *= periPP( 1, 2 ); break; // ep case
case KinematicsMode::ElasticElastic: jacobian_ *= periPP( 2, 2 ); break; // elastic case
case KinematicsMode::InelasticElastic: jacobian_ *= periPP( 3, 2 )*( masses_.dw31_*masses_.dw31_ ); break;
case KinematicsMode::ElasticInelastic: jacobian_ *= periPP( 3, 2 )*( masses_.dw52_*masses_.dw52_ ); break; // single-dissociative case
case KinematicsMode::InelasticInelastic: jacobian_ *= periPP( 3, 3 )*( masses_.dw31_*masses_.dw31_ )*( masses_.dw52_*masses_.dw52_ ); break; // double-dissociative case
}
//--- compute the event weight using the Jacobian
return Constants::GeV2toBarn*jacobian_;
}
//---------------------------------------------------------------------------------------------
void
GamGamLL::fillKinematics( bool )
{
const Particle::Momentum cm = event_->getOneByRole( Particle::IncomingBeam1 ).momentum() + event_->getOneByRole( Particle::IncomingBeam2 ).momentum();
const double gamma = cm.energy()/sqs_, betgam = cm.pz()/sqs_;
Particle::Momentum plab_ip1 = Particle::Momentum( 0., 0., p_cm_, ep1_ ).betaGammaBoost( gamma, betgam );
Particle::Momentum plab_ip2 = Particle::Momentum( 0., 0., -p_cm_, ep2_ ).betaGammaBoost( gamma, betgam );
p3_lab_.betaGammaBoost( gamma, betgam );
p5_lab_.betaGammaBoost( gamma, betgam );
//----- parameterise a random rotation around z-axis
const int rany = ( rand() >= 0.5 * RAND_MAX ) ? 1 : -1,
ransign = ( rand() >= 0.5 * RAND_MAX ) ? 1 : -1;
const double ranphi = ( 0.5 * rand()/RAND_MAX )*2.*M_PI;
Particle::Momentum plab_ph1 = ( plab_ip1-p3_lab_ ).rotatePhi( ranphi, rany );
Particle::Momentum plab_ph2 = ( plab_ip2-p5_lab_ ).rotatePhi( ranphi, rany );
p3_lab_.rotatePhi( ranphi, rany );
p5_lab_.rotatePhi( ranphi, rany );
p6_cm_.rotatePhi( ranphi, rany );
p7_cm_.rotatePhi( ranphi, rany );
/*if ( symmetrise_ && rand() >= .5*RAND_MAX ) {
p6_cm_.mirrorZ();
p7_cm_.mirrorZ();
}*/
//----- incoming protons
event_->getOneByRole( Particle::IncomingBeam1 ).setMomentum( plab_ip1 );
event_->getOneByRole( Particle::IncomingBeam2 ).setMomentum( plab_ip2 );
//----- first outgoing proton
Particle& op1 = event_->getOneByRole( Particle::OutgoingBeam1 );
op1.setMomentum( p3_lab_ );
switch ( cuts_.mode ) {
case KinematicsMode::ElasticElastic:
case KinematicsMode::ElasticInelastic:
default:
op1.setStatus( Particle::Status::FinalState ); // stable proton
break;
case KinematicsMode::InelasticElastic:
case KinematicsMode::InelasticInelastic:
op1.setStatus( Particle::Status::Unfragmented ); // fragmenting remnants
op1.setMass( MX_ );
break;
}
//----- second outgoing proton
Particle& op2 = event_->getOneByRole( Particle::OutgoingBeam2 );
op2.setMomentum( p5_lab_ );
switch ( cuts_.mode ) {
case KinematicsMode::ElasticElastic:
case KinematicsMode::InelasticElastic:
default:
op2.setStatus( Particle::Status::FinalState ); // stable proton
break;
case KinematicsMode::ElasticInelastic:
case KinematicsMode::InelasticInelastic:
op2.setStatus( Particle::Status::Unfragmented ); // fragmenting remnants
op2.setMass( MY_ );
break;
}
//----- first incoming photon
Particle& ph1 = event_->getOneByRole( Particle::Parton1 );
ph1.setMomentum( plab_ph1 );
//----- second incoming photon
Particle& ph2 = event_->getOneByRole( Particle::Parton2 );
ph2.setMomentum( plab_ph2 );
Particles& central_system = event_->getByRole( Particle::CentralSystem );
//----- first outgoing lepton
Particle& ol1 = central_system[0];
ol1.setPdgId( ol1.pdgId(), ransign );
ol1.setMomentum( p6_cm_ );
ol1.setStatus( Particle::Status::FinalState );
//----- second outgoing lepton
Particle& ol2 = central_system[1];
ol2.setPdgId( ol2.pdgId(), -ransign );
ol2.setMomentum( p7_cm_ );
ol2.setStatus( Particle::Status::FinalState );
//----- intermediate two-lepton system
event_->getOneByRole( Particle::Intermediate ).setMomentum( p6_cm_+p7_cm_ );
}
//---------------------------------------------------------------------------------------------
double
GamGamLL::periPP( int nup_, int ndown_ )
{
CG_DEBUG_LOOP( "GamGamLL" )
<< " Nup = " << nup_ << "\n\t"
<< "Ndown = " << ndown_;
FormFactors fp1, fp2;
formFactors( -t1_, -t2_, fp1, fp2 );
CG_DEBUG_LOOP( "GamGamLL" )
<< "u1 = " << fp1.FM << "\n\t"
<< "u2 = " << fp1.FE << "\n\t"
<< "v1 = " << fp2.FM << "\n\t"
<< "v2 = " << fp2.FE;
const double qqq = q1dq_*q1dq_,
qdq = 4.*masses_.Ml2_-w4_;
const double t11 = 64. *( bb_*( qqq-g4_-qdq*( t1_+t2_+2.*masses_.Ml2_ ) )-2.*( t1_+2.*masses_.Ml2_ )*( t2_+2.*masses_.Ml2_ )*qqq ) * t1_*t2_, // magnetic-magnetic
t12 = 128.*( -bb_*( dd2_+g6_ )-2.*( t1_+2.*masses_.Ml2_ )*( sa2_*qqq+a6_*a6_ ) ) * t1_, // electric-magnetic
t21 = 128.*( -bb_*( dd4_+g5_ )-2.*( t2_+2.*masses_.Ml2_ )*( sa1_*qqq+a5_*a5_ ) ) * t2_, // magnetic-electric
t22 = 512.*( bb_*( delta_*delta_-gram_ )-pow( epsi_-delta_*( qdq+q1dq2_ ), 2 )-sa1_*a6_*a6_-sa2_*a5_*a5_-sa1_*sa2_*qqq ); // electric-electric
const double peripp = ( fp1.FM*fp2.FM*t11 + fp1.FE*fp2.FM*t21 + fp1.FM*fp2.FE*t12 + fp1.FE*fp2.FE*t22 ) / pow( 2.*t1_*t2_*bb_, 2 );
CG_DEBUG_LOOP( "GamGamLL" )
<< "t11 = " << t11 << "\t"
<< "t12 = " << t12 << "\n\t"
<< "t21 = " << t21 << "\t"
<< "t22 = " << t22 << "\n\t"
<< "=> PeriPP = " << peripp;
return peripp;
}
void
GamGamLL::map( double expo, const Limits& lim, double& out, double& dout, const std::string& var_name_ )
{
const double y = lim.max()/lim.min();
out = lim.min()*pow( y, expo );
dout = out*log( y );
CG_DEBUG_LOOP( "GamGamLL:map" )
<< "Mapping variable \"" << var_name_ << "\"\n\t"
<< "limits = " << lim << "\n\t"
<< "max/min = " << y << "\n\t"
<< "exponent = " << expo << "\n\t"
<< "output = " << out << "\n\t"
<< "d(output) = " << dout;
}
void
GamGamLL::mapla( double y, double z, int u, double xm, double xp, double& x, double& d )
{
const double xmb = xm-y-z, xpb = xp-y-z;
const double c = -4.*y*z;
const double alp = sqrt( xpb*xpb + c ), alm = sqrt( xmb*xmb + c );
const double am = xmb+alm, ap = xpb+alp;
const double yy = ap/am, zz = pow( yy, u );
x = y + z + ( am*zz - c / ( am*zz ) ) / 2.;
const double ax = sqrt( pow( x-y-z, 2 ) + c );
d = ax*log( yy );
}
void
GamGamLL::formFactors( double q1, double q2, FormFactors& fp1, FormFactors& fp2 ) const
{
const double mx2 = MX_*MX_, my2 = MY_*MY_;
switch ( cuts_.mode ) {
case KinematicsMode::ElectronElectron: default: {
fp1 = FormFactors::trivial(); // electron (trivial) form factor
fp2 = FormFactors::trivial(); // electron (trivial) form factor
} break;
case KinematicsMode::ProtonElectron: {
fp1 = FormFactors::protonElastic( -t1_ ); // proton elastic form factor
fp2 = FormFactors::trivial(); // electron (trivial) form factor
} break;
case KinematicsMode::ElectronProton: {
fp1 = FormFactors::trivial(); // electron (trivial) form factor
fp2 = FormFactors::protonElastic( -t2_ ); // proton elastic form factor
} break;
case KinematicsMode::ElasticElastic: {
fp1 = FormFactors::protonElastic( -t1_ ); // proton elastic form factor
fp2 = FormFactors::protonElastic( -t2_ ); // proton elastic form factor
} break;
case KinematicsMode::ElasticInelastic: {
fp1 = FormFactors::protonElastic( -t1_ );
fp2 = FormFactors::protonInelastic( -t2_, w2_, my2, *cuts_.structure_functions );
} break;
case KinematicsMode::InelasticElastic: {
fp1 = FormFactors::protonInelastic( -t1_, w1_, mx2, *cuts_.structure_functions );
fp2 = FormFactors::protonElastic( -t2_ );
} break;
case KinematicsMode::InelasticInelastic: {
fp1 = FormFactors::protonInelastic( -t1_, w1_, mx2, *cuts_.structure_functions );
fp2 = FormFactors::protonInelastic( -t2_, w2_, my2, *cuts_.structure_functions );
} break;
}
}
// register process and define aliases
REGISTER_PROCESS( lpair, GamGamLL )
REGISTER_PROCESS( gamgamll, GamGamLL )
}
}
-
diff --git a/CepGen/Processes/GamGamLL.h b/CepGen/Processes/GamGamLL.h
index b2039a3..20d7d22 100644
--- a/CepGen/Processes/GamGamLL.h
+++ b/CepGen/Processes/GamGamLL.h
@@ -1,221 +1,221 @@
#ifndef CepGen_Processes_GamGamLL_h
#define CepGen_Processes_GamGamLL_h
#include "CepGen/Processes/GenericProcess.h"
#include "CepGen/Core/ParametersList.h"
namespace CepGen
{
- namespace Process
+ namespace process
{
/**
* Full class of methods and objects to compute the full analytic matrix element
* \cite Vermaseren:1982cz for the \f$\gamma\gamma\to\ell^{+}\ell^{-}\f$ process
* according to a set of kinematic constraints provided for the incoming and
* outgoing particles (the Kinematics object).
* The \a f function created by this Process child has its \a _ndim -dimensional
* coordinates mapped as :
* - 0 = \f$t_1\f$, first incoming photon's virtuality
* - 1 = \f$t_2\f$, second incoming photon's virtuality
* - 2 = \f$s_2\f$ mapping
* - 3 = yy4 = \f$\cos\left(\pi x_3\right)\f$ definition
* - 4 = \f$w_4\f$, the two-photon system's invariant mass
* - 5 = xx6 = \f$\frac{1}{2}\left(1-\cos\theta^{\rm CM}_6\right)\f$ definition (3D rotation of the first outgoing lepton with respect to the two-photon centre-of-mass system). If the \a nm_ optimisation flag is set this angle coefficient value becomes
* \f[\frac{1}{2}\left(\frac{a_{\rm map}}{b_{\rm map}}\frac{\beta-1}{\beta+1}+1\right)\f]
* with \f$a_{\rm map}=\frac{1}{2}\left(w_4-t_1-t_2\right)\f$, \f$b_{\rm map}=\frac{1}{2}\sqrt{\left(\left(w_4-t_1-t_2\right)^2-4t_1t_2\right)\left(1-4\frac{w_6}{w_4}\right)}\f$, and \f$\beta=\left(\frac{a_{\rm map}+b_{\rm map}}{a_{\rm map}-b_{\rm map}}\right)^{2x_5-1}\f$
* and the Jacobian element is scaled by a factor \f$\frac{1}{2}\frac{\left(a_{\rm map}^2-b_{\rm map}^2\cos^2\theta^{\rm CM}_6\right)}{a_{\rm map}b_{\rm map}}\log\left(\frac{a_{\rm map}+b_{\rm map}}{a_{\rm map}-b_{\rm map}}\right)\f$
* - 6 = _phicm6_, or \f$\phi_6^{\rm CM}\f$ the rotation angle of the dilepton system in the centre-of-mass
* system
* - 7 = \f$x_q\f$, \f$w_X\f$ mappings, as used in the single- and double-dissociative
* cases only
* \brief Compute the matrix element for a CE \f$\gamma\gamma\to\ell^{+}\ell^{-}\f$
* process
*/
class GamGamLL : public GenericProcess
{
public:
/// \brief Class constructor: set the mandatory parameters before integration and events generation
/// \param[in] params General process parameters (nopt = Optimisation, legacy from LPAIR)
explicit GamGamLL( const ParametersList& params = ParametersList() );
ProcessPtr clone( const ParametersList& params ) const override { return ProcessPtr( new GamGamLL( params ) ); }
void addEventContent() override;
void beforeComputeWeight() override;
double computeWeight() override;
unsigned int numDimensions() const override;
void setKinematics( const Kinematics& cuts ) override;
void fillKinematics( bool ) override;
/// Compute the ougoing proton remnant mass
/// \param[in] x A random number (between 0 and 1)
/// \param[in] outmass The maximal outgoing particles' invariant mass
/// \param[in] lepmass The outgoing leptons' mass
/// \param[out] dw The size of the integration bin
/// \return Mass of the outgoing proton remnant
double computeOutgoingPrimaryParticlesMasses( double x, double outmass, double lepmass, double& dw );
/// Set all the kinematic variables for the outgoing proton remnants, and prepare the hadronisation
/// \param[in] part Particle to "prepare" for the hadronisation to be performed
void prepareHadronisation( Particle *part );
private:
/**
* Calculate energies and momenta of the
* 1st, 2nd (resp. the "proton-like" and the "electron-like" incoming particles),
* 3rd (the "proton-like" outgoing particle),
* 4th (the two-photons central system), and
* 5th (the "electron-like" outgoing particle) particles in the overall centre-of-mass frame.
* \brief Energies/momenta computation for the various particles, in the CM system
* \return Success state of the operation
*/
bool orient();
/**
* Compute the expression of the matrix element squared for the \f$\gamma\gamma\rightarrow\ell^{+}\ell^{-}\f$ process.
* It returns the value of the convolution of the form factor or structure functions with the central two-photons matrix element squared.
* \brief Computes the matrix element squared for the requested process
* \return Full matrix element for the two-photon production of a pair of spin\f$-\frac{1}{2}-\f$point particles.
* It is noted as \f[
* M = \frac{1}{4bt_1 t_2}\sum_{i=1}^2\sum_{j=1}^2 u_i v_j t_{ij} = \frac{1}{4}\frac{u_1 v_1 t_{11}+u_2 v_1 t_{21}+u_1 v_2 t_{12}+u_2 v_2 t_{22}}{t_1 t_2 b}
* \f] where \f$b\f$ = \a bb_ is defined in \a ComputeWeight as : \f[
* b = t_1 t_2+\left(w_{\gamma\gamma}\sin^2{\theta^{\rm CM}_6}+4m_\ell\cos^2{\theta^{\rm CM}_6}\right) p_g^2
* \f]
*/
double periPP( int, int );
/**
* Describe the kinematics of the process \f$p_1+p_2\to p_3+p_4+p_5\f$ in terms of Lorentz-invariant variables.
* These variables (along with others) will then be fed into the \a PeriPP method (thus are essential for the evaluation of the full matrix element).
* \return Success state of the operation
*/
bool pickin();
/// Internal switch for the optimised code version (LPAIR legacy ; unimplemented here)
int n_opt_;
int pair_;
Limits w_limits_;
Limits q2_limits_;
Limits mx_limits_;
struct Masses
{
Masses();
/// squared mass of the first proton-like outgoing particle
double MX2_;
/// squared mass of the second proton-like outgoing particle
double MY2_;
/// squared mass of the outgoing leptons
double Ml2_;
/// \f$\delta_2=m_1^2-m_2^2\f$ as defined in Vermaseren's paper
/// \cite Vermaseren:1982cz for the full definition of this quantity
double w12_;
/// \f$\delta_1=m_3^2-m_1^2\f$ as defined in Vermaseren's paper
/// \cite Vermaseren:1982cz for the full definition of this quantity
double w31_;
double dw31_;
/// \f$\delta_4=m_5^2-m_2^2\f$ as defined in Vermaseren's paper
/// \cite Vermaseren:1982cz for the full definition of this quantity
double w52_;
double dw52_;
};
Masses masses_;
/// energy of the first proton-like incoming particle
double ep1_;
/// energy of the second proton-like incoming particle
double ep2_;
double p_cm_;
/// energy of the two-photon central system
double ec4_;
/// 3-momentum norm of the two-photon central system
double pc4_;
/// mass of the two-photon central system
double mc4_;
/// squared mass of the two-photon central system
double w4_;
/// \f$p_{12} = \frac{1}{2}\left(s-m_{p_1}^2-m_{p_2}^2\right)\f$
double p12_;
double p1k2_, p2k1_;
/// \f$p_{13} = -\frac{1}{2}\left(t_1-m_{p_1}^2-m_{p_3}^2\right)\f$
double p13_;
double p14_, p25_;
double q1dq_, q1dq2_;
double s1_, s2_;
double epsi_;
double g5_, g6_;
double a5_, a6_;
double bb_;
double gram_;
double dd1_, dd2_, dd3_;
/// \f$\delta_5=m_4^2-t_1\f$ as defined in Vermaseren's paper
/// \cite Vermaseren:1982cz for the full definition of this quantity
double dd4_;
double dd5_;
/**
* Invariant used to tame divergences in the matrix element computation. It is defined as
* \f[\Delta = \left(p_1\cdot p_2\right)\left(q_1\cdot q_2\right)-\left(p_1\cdot q_2\right)\left(p_2\cdot q_1\right)\f]
* with \f$p_i, q_i\f$ the 4-momenta associated to the incoming proton-like particle and to the photon emitted from it.
*/
double delta_;
double g4_;
double sa1_, sa2_;
double sl1_;
/// cosine of the polar angle for the two-photons centre-of-mass system
double cos_theta4_;
/// sine of the polar angle for the two-photons centre-of-mass system
double sin_theta4_;
double al4_;
double be4_;
double de3_, de5_;
double pt4_;
/// Kinematics of the first incoming proton
Particle::Momentum p1_lab_;
/// Kinematics of the second incoming proton
Particle::Momentum p2_lab_;
/// Kinematics of the first outgoing proton
Particle::Momentum p3_lab_;
/// Kinematics of the two-photon system (in the two-proton CM)
Particle::Momentum p4_lab_;
/// Kinematics of the second outgoing proton
Particle::Momentum p5_lab_;
/// Kinematics of the first outgoing lepton (in the two-proton CM)
Particle::Momentum p6_cm_;
/// Kinematics of the second outgoing lepton (in the two-proton CM)
Particle::Momentum p7_cm_;
double jacobian_;
private:
/**
* Define modified variables of integration to avoid peaks integrations (see \cite Vermaseren:1982cz 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] lim Min/maximal value of the variable
* \param[out] out The new variable definition
* \param[out] dout The bin width 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, const Limits& lim, 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 );
/// Compute the electric/magnetic form factors for the two considered \f$Q^{2}\f$ momenta transfers
void formFactors( double q1, double q2, FormFactors& fp1, FormFactors& fp2 ) const;
};
}
}
#endif
diff --git a/CepGen/Processes/GenericKTProcess.cpp b/CepGen/Processes/GenericKTProcess.cpp
index df22774..90fb99a 100644
--- a/CepGen/Processes/GenericKTProcess.cpp
+++ b/CepGen/Processes/GenericKTProcess.cpp
@@ -1,370 +1,370 @@
#include "CepGen/Processes/GenericKTProcess.h"
#include "CepGen/Core/Exception.h"
#include "CepGen/Core/ParametersList.h"
#include "CepGen/Event/Event.h"
#include "CepGen/StructureFunctions/StructureFunctions.h"
#include "CepGen/StructureFunctions/SigmaRatio.h"
#include "CepGen/Physics/Constants.h"
#include "CepGen/Physics/KTFlux.h"
#include "CepGen/Physics/FormFactors.h"
#include "CepGen/Physics/PDG.h"
#include <iomanip>
namespace CepGen
{
- namespace Process
+ namespace process
{
GenericKTProcess::GenericKTProcess( const ParametersList& params,
const std::string& name,
const std::string& description,
const std::array<PDG,2>& partons,
const std::vector<PDG>& central ) :
GenericProcess( name, description+" (kT-factorisation approach)" ),
num_dimensions_( 0 ), kt_jacobian_( 0. ),
qt1_( 0. ), phi_qt1_( 0. ), qt2_( 0. ), phi_qt2_( 0. ),
kIntermediateParts( partons ), kProducedParts( central )
{}
void
GenericKTProcess::addEventContent()
{
GenericProcess::setEventContent(
{ // incoming state
{ Particle::IncomingBeam1, PDG::proton },
{ Particle::IncomingBeam2, PDG::proton },
{ Particle::Parton1, kIntermediateParts[0] },
{ Particle::Parton2, kIntermediateParts[1] }
},
{ // outgoing state
{ Particle::OutgoingBeam1, { PDG::proton } },
{ Particle::OutgoingBeam2, { PDG::proton } },
{ Particle::CentralSystem, kProducedParts }
}
);
setExtraContent();
}
unsigned int
GenericKTProcess::numDimensions() const
{
return num_dimensions_;
}
void
GenericKTProcess::setKinematics( const Kinematics& kin )
{
cuts_ = kin;
const KTFlux flux1 = (KTFlux)cuts_.incoming_beams.first.kt_flux,
flux2 = (KTFlux)cuts_.incoming_beams.second.kt_flux;
if ( cuts_.mode == KinematicsMode::invalid ) {
bool el1 = ( flux1 == KTFlux::P_Photon_Elastic
|| flux1 == KTFlux::HI_Photon_Elastic
|| flux1 == KTFlux::P_Gluon_KMR );
bool el2 = ( flux2 == KTFlux::P_Photon_Elastic
|| flux2 == KTFlux::HI_Photon_Elastic
|| flux2 == KTFlux::P_Gluon_KMR );
if ( el1 && el2 )
cuts_.mode = KinematicsMode::ElasticElastic;
else if ( el1 )
cuts_.mode = KinematicsMode::ElasticInelastic;
else if ( el2 )
cuts_.mode = KinematicsMode::InelasticElastic;
else
cuts_.mode = KinematicsMode::InelasticInelastic;
}
else {
//==========================================================================================
// ensure the first incoming flux is compatible with the kinematics mode
//==========================================================================================
if ( ( cuts_.mode == KinematicsMode::ElasticElastic
|| cuts_.mode == KinematicsMode::ElasticInelastic )
&& ( flux1 != KTFlux::P_Photon_Elastic ) ) {
cuts_.incoming_beams.first.kt_flux = ( (HeavyIon)cuts_.incoming_beams.first.pdg )
? KTFlux::HI_Photon_Elastic
: KTFlux::P_Photon_Elastic;
CG_DEBUG( "GenericKTProcess:kinematics" )
<< "Set the kt flux for first incoming photon to \""
<< cuts_.incoming_beams.first.kt_flux << "\".";
}
else if ( flux1 != KTFlux::P_Photon_Inelastic
&& flux1 != KTFlux::P_Photon_Inelastic_Budnev ) {
if ( (HeavyIon)cuts_.incoming_beams.first.pdg )
throw CG_FATAL( "GenericKTProcess:kinematics" )
<< "Inelastic photon emission from HI not yet supported!";
cuts_.incoming_beams.first.kt_flux = KTFlux::P_Photon_Inelastic_Budnev;
CG_DEBUG( "GenericKTProcess:kinematics" )
<< "Set the kt flux for first incoming photon to \""
<< cuts_.incoming_beams.first.kt_flux << "\".";
}
//==========================================================================================
// ensure the second incoming flux is compatible with the kinematics mode
//==========================================================================================
if ( ( cuts_.mode == KinematicsMode::ElasticElastic
|| cuts_.mode == KinematicsMode::InelasticElastic )
&& ( flux2 != KTFlux::P_Photon_Elastic ) ) {
cuts_.incoming_beams.second.kt_flux = ( (HeavyIon)cuts_.incoming_beams.second.pdg )
? KTFlux::HI_Photon_Elastic
: KTFlux::P_Photon_Elastic;
CG_DEBUG( "GenericKTProcess:kinematics" )
<< "Set the kt flux for second incoming photon to \""
<< cuts_.incoming_beams.second.kt_flux << "\".";
}
else if ( flux2 != KTFlux::P_Photon_Inelastic
&& flux2 != KTFlux::P_Photon_Inelastic_Budnev ) {
if ( (HeavyIon)cuts_.incoming_beams.second.pdg )
throw CG_FATAL( "GenericKTProcess:kinematics" )
<< "Inelastic photon emission from HI not yet supported!";
cuts_.incoming_beams.second.kt_flux = KTFlux::P_Photon_Inelastic_Budnev;
CG_DEBUG( "GenericKTProcess:kinematics" )
<< "Set the kt flux for second incoming photon to \""
<< cuts_.incoming_beams.second.kt_flux << "\".";
}
}
//============================================================================================
// initialise the "constant" (wrt x) part of the Jacobian
//============================================================================================
kt_jacobian_ = 1.;
num_dimensions_ = 0;
mapped_variables_.clear();
//============================================================================================
// register the incoming partons' variables
//============================================================================================
registerVariable( qt1_, Mapping::logarithmic, cuts_.cuts.initial.qt, { 1.e-10, 500. }, "First incoming parton virtuality" );
registerVariable( qt2_, Mapping::logarithmic, cuts_.cuts.initial.qt, { 1.e-10, 500. }, "Second incoming parton virtuality" );
registerVariable( phi_qt1_, Mapping::linear, cuts_.cuts.initial.phi_qt, { 0., 2.*M_PI }, "First incoming parton azimuthal angle" );
registerVariable( phi_qt2_, Mapping::linear, cuts_.cuts.initial.phi_qt, { 0., 2.*M_PI }, "Second incoming parton azimuthal angle" );
//============================================================================================
// register all process-dependent variables
//============================================================================================
preparePhaseSpace();
//============================================================================================
// register the outgoing remnants' variables
//============================================================================================
MX_ = MY_ = event_->getOneByRole( Particle::IncomingBeam1 ).mass();
if ( cuts_.mode == KinematicsMode::InelasticElastic || cuts_.mode == KinematicsMode::InelasticInelastic )
registerVariable( MX_, Mapping::square, cuts_.cuts.remnants.mass_single, { 1.07, 1000. }, "Positive z proton remnant mass" );
if ( cuts_.mode == KinematicsMode::ElasticInelastic || cuts_.mode == KinematicsMode::InelasticInelastic )
registerVariable( MY_, Mapping::square, cuts_.cuts.remnants.mass_single, { 1.07, 1000. }, "Negative z proton remnant mass" );
prepareKinematics();
}
double
GenericKTProcess::computeWeight()
{
if ( mapped_variables_.size() == 0 )
throw CG_FATAL( "GenericKTProcess:weight" )
<< "No variables are mapped with this process!";
if ( kt_jacobian_ == 0. )
throw CG_FATAL( "GenericKTProcess:weight" )
<< "Point-independant component of the Jacobian for this "
<< "kt-factorised process is null.\n\t"
<< "Please check the validity of the phase space!";
//============================================================================================
// generate and initialise all variables, and auxiliary (x-dependent) part of the Jacobian
// for this phase space point.
//============================================================================================
const double aux_jacobian = generateVariables();
if ( aux_jacobian <= 0. )
return 0.;
//============================================================================================
// compute the integrand and combine together into a single weight for the phase space point.
//============================================================================================
const double integrand = computeKTFactorisedMatrixElement();
if ( integrand <= 0. )
return 0.;
const double weight = ( kt_jacobian_*aux_jacobian ) * integrand;
CG_DEBUG_LOOP( "GenericKTProcess:weight" )
<< "Jacobian: " << kt_jacobian_ << " * " << aux_jacobian
<< " = " << ( kt_jacobian_*aux_jacobian ) << ".\n\t"
<< "Integrand = " << integrand << "\n\t"
<< "dW = " << weight << ".";
return weight;
}
std::pair<double,double>
GenericKTProcess::incomingFluxes( double x1, double q1t2, double x2, double q2t2 ) const
{
//--- compute fluxes according to modelling specified in parameters card
std::pair<double,double> fluxes = {
ktFlux( (KTFlux)cuts_.incoming_beams.first.kt_flux, x1, q1t2, *cuts_.structure_functions, MX_ ),
ktFlux( (KTFlux)cuts_.incoming_beams.second.kt_flux, x2, q2t2, *cuts_.structure_functions, MY_ )
};
CG_DEBUG_LOOP( "GenericKTProcess:fluxes" )
<< "KT fluxes: " << fluxes.first << " / " << fluxes.second << ".";
return fluxes;
}
void
GenericKTProcess::registerVariable( double& out, const Mapping& type,
const Limits& in, Limits default_limits,
const char* description )
{
Limits lim = in;
out = 0.; // reset the variable
if ( !in.valid() ) {
CG_DEBUG( "GenericKTProcess:registerVariable" )
<< description << " could not be retrieved from the user configuration!\n\t"
<< "Setting it to the default value: " << default_limits << ".";
lim = default_limits;
}
if ( type == Mapping::logarithmic )
lim = {
std::max( log( lim.min() ), -10. ),
std::min( log( lim.max() ), +10. )
};
mapped_variables_.emplace_back( MappingVariable{ description, lim, out, type, num_dimensions_++ } );
switch ( type ) {
case Mapping::square:
kt_jacobian_ *= 2.*lim.range();
break;
default:
kt_jacobian_ *= lim.range();
break;
}
CG_DEBUG( "GenericKTProcess:registerVariable" )
<< description << " has been mapped to variable " << num_dimensions_ << ".\n\t"
<< "Allowed range for integration: " << lim << ".\n\t"
<< "Variable integration mode: " << type << ".";
}
void
GenericKTProcess::dumpVariables() const
{
std::ostringstream os;
for ( const auto& var : mapped_variables_ )
os << "\n\t(" << var.index << ") " << var.type << " mapping (" << var.description << ") in range " << var.limits;
CG_INFO( "GenericKTProcess:dumpVariables" )
<< "List of variables handled by this kt-factorised process:"
<< os.str();
}
double
GenericKTProcess::generateVariables() const
{
double jacobian = 1.;
for ( const auto& cut : mapped_variables_ ) {
if ( !cut.limits.valid() )
continue;
const double xv = x( cut.index ); // between 0 and 1
switch ( cut.type ) {
case Mapping::linear: {
cut.variable = cut.limits.x( xv );
} break;
case Mapping::logarithmic: {
cut.variable = exp( cut.limits.x( xv ) );
jacobian *= cut.variable;
} break;
case Mapping::square: {
cut.variable = cut.limits.x( xv );
jacobian *= cut.variable;
} break;
}
}
if ( CG_EXCEPT_MATCH( "KtProcess:vars", debugInsideLoop ) ) {
std::ostringstream oss;
for ( const auto& cut : mapped_variables_ ) {
oss << "variable " << cut.index
<< " in range " << std::left << std::setw( 20 ) << cut.limits << std::right
<< " has value " << cut.variable << "\n\t";
}
CG_DEBUG_LOOP( "KtProcess:vars" ) << oss.str();
}
return jacobian;
}
void
GenericKTProcess::fillKinematics( bool )
{
fillCentralParticlesKinematics(); // process-dependent!
fillPrimaryParticlesKinematics();
}
void
GenericKTProcess::fillPrimaryParticlesKinematics()
{
//============================================================================================
// outgoing protons
//============================================================================================
Particle& op1 = event_->getOneByRole( Particle::OutgoingBeam1 ),
&op2 = event_->getOneByRole( Particle::OutgoingBeam2 );
op1.setMomentum( PX_ );
op2.setMomentum( PY_ );
switch ( cuts_.mode ) {
case KinematicsMode::ElasticElastic:
op1.setStatus( Particle::Status::FinalState );
op2.setStatus( Particle::Status::FinalState );
break;
case KinematicsMode::ElasticInelastic:
op1.setStatus( Particle::Status::FinalState );
op2.setStatus( Particle::Status::Unfragmented ); op2.setMass( MY_ );
break;
case KinematicsMode::InelasticElastic:
op1.setStatus( Particle::Status::Unfragmented ); op1.setMass( MX_ );
op2.setStatus( Particle::Status::FinalState );
break;
case KinematicsMode::InelasticInelastic:
op1.setStatus( Particle::Status::Unfragmented ); op1.setMass( MX_ );
op2.setStatus( Particle::Status::Unfragmented ); op2.setMass( MY_ );
break;
default: {
throw CG_FATAL( "GenericKTProcess" )
<< "This kT factorisation process is intended for p-on-p collisions! Aborting.";
} break;
}
//============================================================================================
// incoming partons (photons, pomerons, ...)
//============================================================================================
Particle& g1 = event_->getOneByRole( Particle::Parton1 );
g1.setMomentum( event_->getOneByRole( Particle::IncomingBeam1 ).momentum()-PX_, true );
Particle& g2 = event_->getOneByRole( Particle::Parton2 );
g2.setMomentum( event_->getOneByRole( Particle::IncomingBeam2 ).momentum()-PY_, true );
//============================================================================================
// two-parton system
//============================================================================================
event_->getOneByRole( Particle::Intermediate ).setMomentum( g1.momentum()+g2.momentum() );
}
std::ostream&
operator<<( std::ostream& os, const GenericKTProcess::Mapping& type )
{
switch ( type ) {
case GenericKTProcess::Mapping::linear: return os << "linear";
case GenericKTProcess::Mapping::logarithmic: return os << "logarithmic";
case GenericKTProcess::Mapping::square: return os << "squared";
}
return os;
}
}
}
diff --git a/CepGen/Processes/GenericKTProcess.h b/CepGen/Processes/GenericKTProcess.h
index a1fd04c..2a99ddd 100644
--- a/CepGen/Processes/GenericKTProcess.h
+++ b/CepGen/Processes/GenericKTProcess.h
@@ -1,143 +1,142 @@
#ifndef CepGen_Processes_GenericKTProcess_h
#define CepGen_Processes_GenericKTProcess_h
#include "GenericProcess.h"
namespace CepGen
{
class ParametersList;
- namespace Process
+ namespace process
{
/**
* A generic kT-factorisation process.
* \note
* - First 4 dimensions of the phase space are required for the
* incoming partons' virtualities (radial and azimuthal coordinates).
* - Last 0-2 dimensions may be used for the scattered diffractive
* system(s)' invariant mass definition.
* \brief Class template to define any kT-factorisation process
* \author Laurent Forthomme <laurent.forthomme@cern.ch>
* \date Apr 2016
*/
class GenericKTProcess : public GenericProcess
{
public:
/// Class constructor
/// \param[in] params Parameters list
/// \param[in] name Generic process name
/// \param[in] description Human-readable kT-factorised process name
/// \param[in] partons First and second incoming parton
/// \param[in] output Produced final state particles
GenericKTProcess( const ParametersList& params,
const std::string& name,
const std::string& description,
const std::array<PDG,2>& partons,
const std::vector<PDG>& output );
/// Populate the event content with the generated process' topology
void addEventContent() override;
/// Retrieve the total number of dimensions on which the integration is being performet
unsigned int numDimensions() const override;
/// Retrieve the event weight in the phase space
double computeWeight() override;
/// Populate the event content with the generated process' kinematics
void fillKinematics( bool ) override;
/// List all variables handled by this generic process
void dumpVariables() const;
protected:
/// Set the kinematics associated to the phase space definition
void setKinematics( const Kinematics& kin ) override;
/// Set the kinematics of the central system before any point computation
virtual void setExtraContent() {}
/// Prepare the central part of the Jacobian (only done once, as soon as the kinematics is set)
virtual void preparePhaseSpace() = 0;
/// kT-factorised matrix element (event weight)
/// \return Weight of the point in the phase space to the integral
virtual double computeKTFactorisedMatrixElement() = 0;
/// Compute the unintegrated photon fluxes (for inelastic distributions, interpolation on double logarithmic grid)
std::pair<double,double> incomingFluxes( double, double, double, double ) const;
/// Set the kinematics of the incoming and outgoing protons (or remnants)
void fillPrimaryParticlesKinematics();
/// Set the kinematics of the outgoing central system
virtual void fillCentralParticlesKinematics() = 0;
/// Type of mapping to apply on the variable
enum class Mapping
{
/// a linear \f${\rm d}x\f$ mapping
linear = 0,
/// a logarithmic \f$\frac{{\rm d}x}{x} = {\rm d}(\log x)\f$ mapping
logarithmic,
/// a square \f${\rm d}x^2=2x\cdot{\rm d}x\f$ mapping
square
};
friend std::ostream& operator<<( std::ostream&, const Mapping& );
/// Register a variable to be handled and populated whenever
/// a new phase space point weight is to be calculated.
/// \note To be run once per generation (before any point computation)
/// \param[out] out Reference to the variable to be mapped
/// \param[in] type Type of mapping to apply
/// \param[in] in Integration limits
/// \param[in] default_limits Limits to apply if none retrieved from the user configuration
/// \param[in] description Human-readable description of the variable
void registerVariable( double& out, const Mapping& type, const Limits& in, Limits default_limits, const char* description );
/// Generate and initialise all variables handled by this process
/// \return Phase space point-dependent component of the Jacobian weight of the point in the phase space for integration
/// \note To be run at each point computation (therefore, to be optimised!)
double generateVariables() const;
/// Number of dimensions on which to perform the integration
unsigned short num_dimensions_;
/// Phase space point-independant component of the Jacobian weight of the point in the phase space for integration
double kt_jacobian_;
/// Log-virtuality range of the intermediate parton
Limits log_qt_limits_;
/// Intermediate azimuthal angle range
Limits phi_qt_limits_;
/// Invariant mass range for the scattered excited system
Limits mx_limits_;
/// Virtuality of the first intermediate parton (photon, pomeron, ...)
double qt1_;
/// Azimuthal rotation of the first intermediate parton's transverse virtuality
double phi_qt1_;
/// Virtuality of the second intermediate parton (photon, pomeron, ...)
double qt2_;
/// Azimuthal rotation of the second intermediate parton's transverse virtuality
double phi_qt2_;
/// First outgoing proton
Particle::Momentum PX_;
/// Second outgoing proton
Particle::Momentum PY_;
/// Handler to a variable mapped by this process
struct MappingVariable
{
/// Human-readable description of the variable
std::string description;
/// Kinematic limits to apply on the variable
Limits limits;
/// Reference to the process variable to generate/map
double& variable;
/// Interpolation type
Mapping type;
/// Corresponding integration variable
unsigned short index;
};
/// Collection of variables to be mapped at the weight generation stage
std::vector<MappingVariable> mapped_variables_;
private:
/// First and second intermediate parton (photon, pomeron, ...)
std::array<PDG,2> kIntermediateParts;
/// Type of particles produced in the final state
std::vector<PDG> kProducedParts;
};
}
}
#endif
-
diff --git a/CepGen/Processes/GenericProcess.cpp b/CepGen/Processes/GenericProcess.cpp
index 84fe7c4..7f3f719 100644
--- a/CepGen/Processes/GenericProcess.cpp
+++ b/CepGen/Processes/GenericProcess.cpp
@@ -1,228 +1,228 @@
#include "CepGen/Processes/GenericProcess.h"
#include "CepGen/Core/Exception.h"
#include "CepGen/Event/Event.h"
#include "CepGen/Physics/ParticleProperties.h"
#include "CepGen/Physics/Constants.h"
#include "CepGen/Physics/FormFactors.h"
#include "CepGen/Physics/PDG.h"
namespace CepGen
{
- namespace Process
+ namespace process
{
const double GenericProcess::mp_ = ParticleProperties::mass( PDG::proton );
const double GenericProcess::mp2_ = GenericProcess::mp_*GenericProcess::mp_;
GenericProcess::GenericProcess( const std::string& name, const std::string& description, bool has_event ) :
name_( name ), description_( description ),
first_run( true ),
s_( -1. ), sqs_( -1. ),
MX_( -1. ), MY_( -1. ), w1_( -1. ), w2_( -1. ),
t1_( -1. ), t2_( -1. ),
has_event_( has_event ), event_( new Event ),
is_point_set_( false )
{}
GenericProcess::GenericProcess( const GenericProcess& proc ) :
name_( proc.name_ ), description_( proc.description_ ),
first_run( proc.first_run ),
s_( proc.s_ ), sqs_( proc.sqs_ ),
MX_( proc.MX_ ), MY_( proc.MY_ ), w1_( proc.w1_ ), w2_( proc.w2_ ),
t1_( -1. ), t2_( -1. ), cuts_( proc.cuts_ ),
has_event_( proc.has_event_ ), event_( new Event( *proc.event_.get() ) ),
is_point_set_( false )
{}
GenericProcess&
GenericProcess::operator=( const GenericProcess& proc )
{
name_ = proc.name_; description_ = proc.description_;
first_run = proc.first_run;
s_ = proc.s_; sqs_ = proc.sqs_;
MX_ = proc.MX_; MY_ = proc.MY_; w1_ = proc.w1_; w2_ = proc.w2_;
cuts_ = proc.cuts_;
has_event_ = proc.has_event_; event_.reset( new Event( *proc.event_.get() ) );
is_point_set_ = false;
return *this;
}
void
GenericProcess::setPoint( const unsigned int ndim, double* x )
{
x_ = std::vector<double>( x, x+ndim );
is_point_set_ = true;
if ( CG_EXCEPT_MATCH( "Process:dumpPoint", debugInsideLoop ) )
dumpPoint();
}
double
GenericProcess::x( unsigned int idx ) const
{
if ( idx >= x_.size() )
return -1.;
return x_[idx];
}
void
GenericProcess::clearEvent()
{
event_->restore();
}
void
GenericProcess::setKinematics( const Kinematics& cuts )
{
cuts_ = cuts;
prepareKinematics();
}
void
GenericProcess::prepareKinematics()
{
if ( !isKinematicsDefined() )
throw CG_FATAL( "GenericProcess" ) << "Kinematics not properly defined for the process.";
// at some point introduce non head-on colliding beams?
Particle::Momentum p1( 0., 0., cuts_.incoming_beams.first.pz ), p2( 0., 0., -cuts_.incoming_beams.second.pz );
// on-shell beam particles
p1.setMass( ParticleProperties::mass( cuts_.incoming_beams.first.pdg ) );
p2.setMass( ParticleProperties::mass( cuts_.incoming_beams.second.pdg ) );
setIncomingKinematics( p1, p2 );
sqs_ = CMEnergy( p1, p2 );
s_ = sqs_*sqs_;
w1_ = p1.mass2();
w2_ = p2.mass2();
CG_DEBUG( "GenericProcess" ) << "Kinematics successfully prepared! sqrt(s) = " << sqs_ << ".";
}
void
GenericProcess::dumpPoint() const
{
std::ostringstream os;
for ( unsigned short i = 0; i < x_.size(); ++i ) {
os << Form( " x(%2d) = %8.6f\n\t", i, x_[i] );
}
CG_INFO( "GenericProcess" )
<< "Number of integration parameters: " << x_.size() << "\n\t"
<< os.str() << ".";
}
void
GenericProcess::setEventContent( const IncomingState& ini, const OutgoingState& fin )
{
if ( !has_event_ )
return;
event_->clear();
//----- add the particles in the event
//--- incoming state
for ( const auto& ip : ini ) {
Particle& p = event_->addParticle( ip.first );
p.setPdgId( ip.second, ParticleProperties::charge( ip.second ) );
p.setMass( ParticleProperties::mass( ip.second ) );
if ( ip.first == Particle::IncomingBeam1
|| ip.first == Particle::IncomingBeam2 )
p.setStatus( Particle::Status::PrimordialIncoming );
if ( ip.first == Particle::Parton1
|| ip.first == Particle::Parton2 )
p.setStatus( Particle::Status::Incoming );
}
//--- central system (if not already there)
const auto& central_system = ini.find( Particle::CentralSystem );
if ( central_system == ini.end() ) {
Particle& p = event_->addParticle( Particle::Intermediate );
p.setPdgId( PDG::invalid );
p.setStatus( Particle::Status::Propagator );
}
//--- outgoing state
for ( const auto& opl : fin ) { // pair(role, list of PDGids)
for ( const auto& pdg : opl.second ) {
Particle& p = event_->addParticle( opl.first );
p.setPdgId( pdg, ParticleProperties::charge( pdg ) );
p.setMass( ParticleProperties::mass( pdg ) );
}
}
//----- define the particles parentage
const Particles parts = event_->particles();
for ( const auto& p : parts ) {
Particle& part = event_->operator[]( p.id() );
switch ( part.role() ) {
case Particle::OutgoingBeam1:
case Particle::Parton1:
part.addMother( event_->getOneByRole( Particle::IncomingBeam1 ) );
break;
case Particle::OutgoingBeam2:
case Particle::Parton2:
part.addMother( event_->getOneByRole( Particle::IncomingBeam2 ) );
break;
case Particle::Intermediate:
part.addMother( event_->getOneByRole( Particle::Parton1 ) );
part.addMother( event_->getOneByRole( Particle::Parton2 ) );
break;
case Particle::CentralSystem:
part.addMother( event_->getOneByRole( Particle::Intermediate ) );
break;
default: break;
}
}
//----- freeze the event as it is
event_->freeze();
last_event = event_;
}
void
GenericProcess::setIncomingKinematics( const Particle::Momentum& p1, const Particle::Momentum& p2 )
{
if ( !has_event_ )
return;
event_->getOneByRole( Particle::IncomingBeam1 ).setMomentum( p1 );
event_->getOneByRole( Particle::IncomingBeam2 ).setMomentum( p2 );
}
bool
GenericProcess::isKinematicsDefined()
{
if ( !has_event_ )
return true;
// check the incoming state
bool is_incoming_state_set =
( !event_->getByRole( Particle::IncomingBeam1 ).empty()
&& !event_->getByRole( Particle::IncomingBeam2 ).empty() );
// check the outgoing state
bool is_outgoing_state_set =
( !event_->getByRole( Particle::OutgoingBeam1 ).empty()
&& !event_->getByRole( Particle::OutgoingBeam2 ).empty()
&& !event_->getByRole( Particle::CentralSystem ).empty() );
// combine both states
return is_incoming_state_set && is_outgoing_state_set;
}
std::ostream&
operator<<( std::ostream& os, const GenericProcess& proc )
{
return os << proc.name().c_str();
}
std::ostream&
operator<<( std::ostream& os, const GenericProcess* proc )
{
return os << proc->name().c_str();
}
}
}
diff --git a/CepGen/Processes/GenericProcess.h b/CepGen/Processes/GenericProcess.h
index 34936b3..135cf14 100644
--- a/CepGen/Processes/GenericProcess.h
+++ b/CepGen/Processes/GenericProcess.h
@@ -1,167 +1,167 @@
#ifndef CepGen_Processes_GenericProcess_h
#define CepGen_Processes_GenericProcess_h
#include "CepGen/Event/Particle.h"
#include "CepGen/Physics/Kinematics.h"
#include <vector>
#include <memory>
namespace CepGen
{
class Event;
class FormFactors;
class ParametersList;
/// Location for all physics processes to be generated
- namespace Process
+ namespace process
{
/// \brief Class template to define any process to compute using this MC integrator/events generator
/// \author Laurent Forthomme <laurent.forthomme@cern.ch>
/// \date Jan 2014
class GenericProcess
{
public:
/// Default constructor for an undefined process
/// \param[in] name Process name
/// \param[in] description Human-readable description of the process
/// \param[in] has_event Do we generate the associated event structure?
GenericProcess( const std::string& name, const std::string& description = "<invalid process>", bool has_event = true );
/// Copy constructor for a user process
GenericProcess( const GenericProcess& );
virtual ~GenericProcess() = default;
/// Assignment operator
GenericProcess& operator=( const GenericProcess& );
/// Human-readable format dump of a GenericProcess object
friend std::ostream& operator<<( std::ostream& os, const GenericProcess& proc );
/// Human-readable format dump of a pointer to a GenericProcess object
friend std::ostream& operator<<( std::ostream& os, const GenericProcess* proc );
/// Generic map of particles with their role in the process
typedef std::map<Particle::Role,PDG> ParticlesRoleMap;
/// Pair of particle with their associated role in the process
typedef std::pair<Particle::Role,PDG> ParticleWithRole;
/// Map of all incoming state particles in the process
typedef ParticlesRoleMap IncomingState;
/// Map of all outgoing particles in the process
typedef std::map<Particle::Role,std::vector<PDG> > OutgoingState;
/// Copy all process attributes into a new object
virtual std::unique_ptr<GenericProcess> clone( const ParametersList& ) const = 0;
/// Restore the Event object to its initial state
void clearEvent();
/// Set the kinematics of the incoming state particles
void setIncomingKinematics( const Particle::Momentum& p1, const Particle::Momentum& p2 );
/// Compute the incoming state kinematics
void prepareKinematics();
public:
/// Set the incoming and outgoing state to be expected in the process
inline virtual void addEventContent() {}
/// Set the list of kinematic cuts to apply on the outgoing particles' final state
/// \param[in] cuts The Cuts object containing the kinematic parameters
virtual void setKinematics( const Kinematics& cuts );
/// Return the number of dimensions on which the integration has to be performed
/// \return Number of dimensions on which to integrate
virtual unsigned int numDimensions() const = 0;
/// Prepare the process for its integration over the whole phase space
inline virtual void beforeComputeWeight() {}
/// Compute the weight for this point in the phase-space
virtual double computeWeight() = 0;
/// Fill the Event object with the particles' kinematics
/// \param[in] symmetrise Symmetrise the event? (randomise the production of positively- and negatively-charged outgoing central particles)
virtual void fillKinematics( bool symmetrise = false ) = 0;
public:
/**
* Sets the phase space point to compute the weight associated to it.
* \brief Sets the phase space point to compute
* \param[in] ndim The number of dimensions of the point in the phase space
* \param[in] x[] The (\a ndim_)-dimensional point in the phase space on which the kinematics and the cross-section are computed
*/
void setPoint( const unsigned int ndim, double* x );
/// Dump the evaluated point's coordinates in the standard output stream
void dumpPoint() const;
/// Complete list of Particle with their role in the process for the point considered in the phase space, returned as an Event object.
/// \return Event object containing all the generated Particle objects
inline std::shared_ptr<Event> event() const { return event_; }
///Get the number of dimensions on which the integration is performed
inline const unsigned int ndim() const { return x_.size(); }
/// Get the value of a component of the d-dimensional point considered
double x( unsigned int idx ) const;
/// Name of the process considered
inline const std::string& name() const { return name_; }
/// Human-readable description of the process
inline const std::string& description() const { return description_; }
/// Does the process contain (and hold) an event?
bool hasEvent() const { return has_event_; }
/// Pointer to the last event produced in this run
std::shared_ptr<Event> last_event;
protected:
static const double mp_; ///< Proton mass, in GeV/c\f${}^2\f$
static const double mp2_; ///< Squared proton mass, in GeV\f${}^2\f$/c\f${}^4\f$
/// Set the incoming and outgoing states to be defined in this process (and prepare the Event object accordingly)
void setEventContent( const IncomingState& ini, const OutgoingState& fin );
// ---
/// Name of the process
std::string name_;
/// Process human-readable description
std::string description_;
public:
/// Is it the first time the process is computed?
bool first_run;
protected:
/// Array of double precision floats representing the point on which the weight in the cross-section is computed
std::vector<double> x_;
/// \f$s\f$, squared centre of mass energy of the incoming particles' system, in \f$\mathrm{GeV}^2\f$
double s_;
/// \f$\sqrt s\f$, centre of mass energy of the incoming particles' system (in GeV)
double sqs_;
/// Invariant mass of the first proton-like outgoing particle (or remnant)
double MX_;
/// Invariant mass of the second proton-like outgoing particle (or remnant)
double MY_;
/// \f$m_1^2\f$, squared mass of the first proton-like incoming particle
double w1_;
/// \f$m_2^2\f$, squared mass of the second proton-like incoming particle
double w2_;
/// Virtuality of the first incoming photon
double t1_;
/// Virtuality of the second incoming photon
double t2_;
/// Set of cuts to apply on the final phase space
Kinematics cuts_;
/// Does the process contain (and hold) an event?
bool has_event_;
/// Event object containing all the information on the in- and outgoing particles
std::shared_ptr<Event> event_;
/// Is the phase space point set?
bool is_point_set_;
private:
/**
* Is the system's kinematics well defined and compatible with the process ?
* This check is mandatory to perform the d-dimensional point's cross-section computation.
* \brief Is the system's kinematics well defined?
* \return A boolean stating if the input kinematics and the final states are well-defined
*/
bool isKinematicsDefined();
};
}
/// Helper typedef for a Process unique pointer
- typedef std::unique_ptr<Process::GenericProcess> ProcessPtr;
+ typedef std::unique_ptr<process::GenericProcess> ProcessPtr;
}
#endif
diff --git a/CepGen/Processes/PPtoFF.cpp b/CepGen/Processes/PPtoFF.cpp
index f876e49..ab21577 100644
--- a/CepGen/Processes/PPtoFF.cpp
+++ b/CepGen/Processes/PPtoFF.cpp
@@ -1,379 +1,378 @@
#include "CepGen/Processes/PPtoFF.h"
#include "CepGen/Event/Event.h"
#include "CepGen/Physics/Constants.h"
#include "CepGen/Physics/FormFactors.h"
#include "CepGen/Physics/PDG.h"
#include "CepGen/Core/Exception.h"
#include <iomanip>
#include "CepGen/Processes/ProcessesHandler.h"
namespace CepGen
{
- namespace Process
+ namespace process
{
PPtoFF::PPtoFF( const ParametersList& params ) :
GenericKTProcess( params, "pptoff", "ɣɣ → f⁺f¯", { { PDG::photon, PDG::photon } }, { PDG::muon, PDG::muon } ),
pair_( params.get<int>( "pair", 0 ) ),
method_( params.get<int>( "method", 1 ) ),
y1_( 0. ), y2_( 0. ), pt_diff_( 0. ), phi_pt_diff_( 0. )
{}
void
PPtoFF::preparePhaseSpace()
{
registerVariable( y1_, Mapping::linear, cuts_.cuts.central.rapidity_single, { -6., 6. }, "First outgoing fermion rapidity" );
registerVariable( y2_, Mapping::linear, cuts_.cuts.central.rapidity_single, { -6., 6. }, "Second outgoing fermion rapidity" );
registerVariable( pt_diff_, Mapping::linear, cuts_.cuts.central.pt_diff, { 0., 50. }, "Fermions transverse momentum difference" );
registerVariable( phi_pt_diff_, Mapping::linear, cuts_.cuts.central.phi_pt_diff, { 0., 2.*M_PI }, "Fermions azimuthal angle difference" );
if ( (PDG)pair_ == PDG::invalid )
throw CG_FATAL( "PPtoFF:prepare" )
<< "Invalid fermion pair selected: " << pair_ << "!";
const PDG pdg_f = (PDG)pair_;
mf_ = ParticleProperties::mass( pdg_f );
mf2_ = mf_*mf_;
qf_ = ParticleProperties::charge( pdg_f );
colf_ = ParticleProperties::colours( pdg_f );
CG_DEBUG( "PPtoFF:prepare" )
<< "Produced particles (" << pdg_f << ") "
<< "with mass = " << mf_ << " GeV, "
<< "and charge = " << std::setprecision( 2 ) << qf_ << " e";
CG_DEBUG( "PPtoFF:mode" )
<< "matrix element computation method: " << method_ << ".";
}
double
PPtoFF::computeKTFactorisedMatrixElement()
{
//=================================================================
// matrix element computation
//=================================================================
//--- central partons
const Particle::Momentum q1t( qt1_*cos( phi_qt1_ ), qt1_*sin( phi_qt1_ ), 0. );
const Particle::Momentum q2t( qt2_*cos( phi_qt2_ ), qt2_*sin( phi_qt2_ ), 0. );
CG_DEBUG_LOOP( "PPtoFF" ) << "q(1/2)t = " << q1t << ", " << q2t << ".";
//--- two-parton system
const Particle::Momentum ptsum = q1t+q2t;
const Particle::Momentum ptdiff( pt_diff_*cos( phi_pt_diff_ ), pt_diff_*sin( phi_pt_diff_ ), 0. );
//--- outgoing fermions
const Particle::Momentum p1_cm = 0.5*( ptsum+ptdiff ), p2_cm = 0.5*( ptsum-ptdiff );
//=================================================================
// a window in single particle transverse momentum
//=================================================================
const Limits& pt_limits = cuts_.cuts.central.pt_single;
if ( !pt_limits.passes( p1_cm.pt() ) || !pt_limits.passes( p2_cm.pt() ) )
return 0.;
//=================================================================
// a window in transverse momentum difference
//=================================================================
if ( !cuts_.cuts.central.pt_diff.passes( fabs( p1_cm.pt()-p2_cm.pt() ) ) )
return 0.;
//=================================================================
// a window in rapidity distance
//=================================================================
if ( !cuts_.cuts.central.rapidity_diff.passes( fabs( y1_-y2_ ) ) )
return 0.;
//=================================================================
// auxiliary quantities
//=================================================================
// transverse mass for the two fermions
const double amt1 = std::hypot( p1_cm.pt(), mf_ ), amt2 = std::hypot( p2_cm.pt(), mf_ );
const double alpha1 = amt1/sqs_*exp( y1_ ), beta1 = amt1/sqs_*exp( -y1_ ),
alpha2 = amt2/sqs_*exp( y2_ ), beta2 = amt2/sqs_*exp( -y2_ );
CG_DEBUG_LOOP( "PPtoFF" )
<< "Sudakov parameters:\n\t"
<< " alpha(1/2) = " << alpha1 << ", " << alpha2 << "\n\t"
<< " beta(1/2) = " << beta1 << ", " << beta2 << ".";
const double x1 = alpha1+alpha2, x2 = beta1+beta2;
if ( x1 <= 0. || x1 > 1. || x2 <= 0. || x2 > 1. )
return 0.; // sanity check
//=================================================================
// additional conditions for energy-momentum conservation
//=================================================================
const double s1_eff = x1*s_-qt1_*qt1_, s2_eff = x2*s_-qt2_*qt2_;
const double invm = sqrt( amt1*amt1 + amt2*amt2 + 2.*amt1*amt2*cosh(y1_-y2_) - ptsum.pt2() );
CG_DEBUG_LOOP( "PPtoFF" )
<< "s(1/2)eff = " << s1_eff << ", " << s2_eff << " GeV²\n\t"
<< "central system's invariant mass = " << invm << " GeV.";
if ( ( cuts_.mode == KinematicsMode::ElasticInelastic
|| cuts_.mode == KinematicsMode::InelasticInelastic )
&& ( sqrt( s1_eff ) <= ( MY_+invm ) ) )
return 0.;
if ( ( cuts_.mode == KinematicsMode::InelasticElastic
|| cuts_.mode == KinematicsMode::InelasticInelastic )
&& ( sqrt( s2_eff ) <= ( MX_+invm ) ) )
return 0.;
//=================================================================
// four-momenta of the outgoing protons (or remnants)
//=================================================================
const Particle::Momentum& ak1 = event_->getOneByRole( Particle::IncomingBeam1 ).momentum(),
&ak2 = event_->getOneByRole( Particle::IncomingBeam2 ).momentum();
CG_DEBUG_LOOP( "PPtoFF" )
<< "incoming particles: p(1/2) = " << ak1 << ", " << ak2 << ".";
const double px_plus = ( 1.-x1 )*M_SQRT2*ak1.p(),
px_minus = ( MX_*MX_ + q1t.pt2() )*0.5/px_plus;
const double py_minus = ( 1.-x2 )*M_SQRT2*ak2.p(),
py_plus = ( MY_*MY_ + q2t.pt2() )*0.5/py_minus;
CG_DEBUG_LOOP( "PPtoFF" )
<< "px± = " << px_plus << ", " << px_minus << "\n\t"
<< "py± = " << py_plus << ", " << py_minus << ".";
PX_ = Particle::Momentum( -q1t.px(), -q1t.py(), ( px_plus-px_minus )*M_SQRT1_2, ( px_plus+px_minus )*M_SQRT1_2 );
PY_ = Particle::Momentum( -q2t.px(), -q2t.py(), ( py_plus-py_minus )*M_SQRT1_2, ( py_plus+py_minus )*M_SQRT1_2 );
CG_DEBUG_LOOP( "PPtoFF" )
<< "First remnant: " << PX_ << ", mass = " << PX_.mass() << "\n\t"
<< "Second remnant: " << PY_ << ", mass = " << PY_.mass() << ".";
if ( fabs( PX_.mass()-MX_ ) > 1.e-6 )
throw CG_FATAL( "PPtoFF" ) << "Invalid X system mass: " << PX_.mass() << "/" << MX_ << ".";
if ( fabs( PY_.mass()-MY_ ) > 1.e-6 )
throw CG_FATAL( "PPtoFF" ) << "Invalid Y system mass: " << PY_.mass() << "/" << MY_ << ".";
//=================================================================
// four-momenta of the outgoing l^+ and l^-
//=================================================================
const Particle::Momentum p1 = p1_cm + alpha1*ak1 + beta1*ak2;
const Particle::Momentum p2 = p2_cm + alpha2*ak1 + beta2*ak2;
CG_DEBUG_LOOP( "PPtoFF" )
<< "unboosted first fermion: " << p1 << ", mass = " << p1.mass() << "\n\t"
<< " second fermion: " << p2 << ", mass = " << p2.mass() << ".";
p_f1_ = Particle::Momentum::fromPxPyYM( p1_cm.px(), p1_cm.py(), y2_, mf_ );
p_f2_ = Particle::Momentum::fromPxPyYM( p2_cm.px(), p2_cm.py(), y1_, mf_ );
CG_DEBUG_LOOP( "PPtoFF" )
<< "First fermion: " << p_f1_ << ", mass = " << p_f1_.mass() << "\n\t"
<< "Second fermion: " << p_f2_ << ", mass = " << p_f2_.mass() << ".";
if ( fabs( p_f1_.mass()-mf_ ) > 1.e-4 )
throw CG_FATAL( "PPtoFF" ) << "Invalid fermion 1 mass: " << p_f1_.mass() << "/" << mf_ << ".";
if ( fabs( p_f2_.mass()-mf_ ) > 1.e-4 )
throw CG_FATAL( "PPtoFF" ) << "Invalid fermion 2 mass: " << p_f2_.mass() << "/" << mf_ << ".";
//=================================================================
// matrix elements
//=================================================================
double amat2 = 0.;
//=================================================================
// How matrix element is calculated
//=================================================================
if ( method_ == 0 ) {
//===============================================================
// Mendelstam variables
//===============================================================
//const double shat = s_*x1*x2; // approximation
const double shat = ( q1t+q2t ).mass2(); // exact formula
//const double mll = sqrt( shat );
const double that1 = ( q1t-p1 ).mass2(), that2 = ( q2t-p2 ).mass2();
const double uhat1 = ( q1t-p2 ).mass2(), uhat2 = ( q2t-p1 ).mass2();
const double that = 0.5*( that1+that2 ), uhat = 0.5*( uhat1+uhat2 );
amat2 = onShellME( shat, that, uhat );
CG_DEBUG_LOOP( "PPtoFF:onShell" )
<< "that(1/2) = " << that1 << " / " << that2 << "\n\t"
<< "uhat(1/2) = " << uhat1 << " / " << uhat2 << "\n\t"
<< "squared matrix element: " << amat2 << ".";
}
else if ( method_ == 1 ) {
const double t1abs = ( q1t.pt2() + x1*( MX_*MX_-mp2_ )+x1*x1*mp2_ )/( 1.-x1 ),
t2abs = ( q2t.pt2() + x2*( MY_*MY_-mp2_ )+x2*x2*mp2_ )/( 1.-x2 );
const double z1p = alpha1/x1, z1m = alpha2/x1,
z2p = beta1 /x2, z2m = beta2 /x2;
CG_DEBUG_LOOP( "PPtoFF:offShell" )
<< "z(1/2)p = " << z1p << ", " << z2p << "\n\t"
<< "z(1/2)m = " << z1m << ", " << z2m << ".";
amat2 = offShellME( t1abs, t2abs, z1m, z1p, z2m, z2p, q1t, q2t ) * pow( x1*x2*s_, 2 );
}
//============================================
// unintegrated photon distributions
//============================================
const std::pair<double,double> fluxes
= GenericKTProcess::incomingFluxes( x1, q1t.pt2(), x2, q2t.pt2() );
CG_DEBUG_LOOP( "PPtoFF" )
<< "Incoming photon fluxes for (x/kt2) = "
<< "(" << x1 << "/" << q1t.pt2() << "), "
<< "(" << x2 << "/" << q2t.pt2() << "):\n\t"
<< fluxes.first << ", " << fluxes.second << ".";
//=================================================================
// 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 g_em = 4.*M_PI*Constants::alphaEM*qf_*qf_;
const double aintegral = amat2 * colf_ * ( g_em*g_em )
* 1. / pow( 4.*M_PI*( x1*x2*s_ ), 2 )
* fluxes.first*M_1_PI * fluxes.second*M_1_PI * 0.25
* Constants::GeV2toBarn;
//=================================================================
return aintegral*qt1_*qt2_*pt_diff_;
//=================================================================
}
void
PPtoFF::fillCentralParticlesKinematics()
{
// randomise the charge of the outgoing fermions
short sign = ( drand() > 0.5 ) ? +1 : -1;
//=================================================================
// first outgoing fermion
//=================================================================
Particle& of1 = event_->getByRole( Particle::CentralSystem )[0];
of1.setPdgId( of1.pdgId(), sign );
of1.setStatus( Particle::Status::FinalState );
of1.setMomentum( p_f1_ );
//=================================================================
// second outgoing fermion
//=================================================================
Particle& of2 = event_->getByRole( Particle::CentralSystem )[1];
of2.setPdgId( of2.pdgId(), -sign );
of2.setStatus( Particle::Status::FinalState );
of2.setMomentum( p_f2_ );
}
double
PPtoFF::onShellME( double shat, double that, double uhat ) const
{
CG_DEBUG_LOOP( "PPtoFF:onShell" )
<< "shat: " << shat << ", that: " << that << ", uhat: " << uhat << ".";
//=================================================================
// on-shell formula for M^2
//=================================================================
const double ml4 = mf2_*mf2_, ml8 = ml4*ml4;
const double term1 = 6. *ml8,
term2 = -3. *ml4 *that*that,
term3 = -14.*ml4 *that*uhat,
term4 = -3. *ml4 *uhat*uhat,
term5 = mf2_*that*that*that,
term6 = 7.* mf2_*that*that*uhat,
term7 = 7.* mf2_*that*uhat*uhat,
term8 = mf2_*uhat*uhat*uhat,
term9 = -that*that*that*uhat,
term10 = -that*uhat*uhat*uhat;
return -2.*( term1+term2+term3+term4+term5
+term6+term7+term8+term9+term10 )/( pow( ( mf2_-that )*( mf2_-uhat ), 2) );
}
double
PPtoFF::offShellME( double t1abs, double t2abs, double z1m, double z1p, double z2m, double z2p, const Particle::Momentum& q1, const Particle::Momentum& q2 ) const
{
//=================================================================
// Wolfgang's formulae
//=================================================================
const double z1 = z1p*z1m, z2 = z2p*z2m;
const double eps12 = mf2_+z1*t1abs, eps22 = mf2_+z2*t2abs;
const Particle::Momentum ak1 = ( z1m*p_f1_-z1p*p_f2_ ), ak2 = ( z2m*p_f1_-z2p*p_f2_ );
const Particle::Momentum ph_p1 = ak1+z1p*q2, ph_m1 = ak1-z1m*q2;
const Particle::Momentum ph_p2 = ak2+z2p*q1, ph_m2 = ak2-z2m*q1;
const Particle::Momentum phi1(
ph_p1.px()/( ph_p1.pt2()+eps12 )-ph_m1.px()/( ph_m1.pt2()+eps12 ),
ph_p1.py()/( ph_p1.pt2()+eps12 )-ph_m1.py()/( ph_m1.pt2()+eps12 ),
0.,
1./( ph_p1.pt2()+eps12 )-1./( ph_m1.pt2()+eps12 )
);
const Particle::Momentum phi2(
ph_p2.px()/( ph_p2.pt2()+eps22 )-ph_m2.px()/( ph_m2.pt2()+eps22 ),
ph_p2.py()/( ph_p2.pt2()+eps22 )-ph_m2.py()/( ph_m2.pt2()+eps22 ),
0.,
1./( ph_p2.pt2()+eps22 )-1./( ph_m2.pt2()+eps22 )
);
const double dot1 = phi1.threeProduct( q1 )/qt1_, cross1 = phi1.crossProduct( q1 )/qt1_;
const double dot2 = phi2.threeProduct( q2 )/qt2_, cross2 = phi2.crossProduct( q2 )/qt2_;
CG_DEBUG_LOOP( "PPtoFF:offShell" )
<< "phi1 = " << phi1 << "\n\t"
<< "phi2 = " << phi2 << "\n\t"
<< "(dot): " << dot1 << " / " << dot2 << "\n\t"
<< "(cross): " << cross1 << " / " << cross2 << ".";
//=================================================================
// six terms in Wolfgang's formula for
// off-shell gamma gamma --> l^+ l^-
//=================================================================
const unsigned short imat1 = 1, imat2 = 1;
const unsigned short itermLL = 1, itermTT = 1, itermLT = 1, itermtt = 1;
const double aux2_1 = itermLL * ( mf2_ + 4.*z1*z1*t1abs ) * phi1.energy2()
+itermTT * ( ( z1p*z1p + z1m*z1m )*( dot1*dot1 + cross1*cross1 ) )
+itermtt * ( cross1*cross1 - dot1*dot1 )
-itermLT * 4.*z1*( z1p-z1m ) * phi1.energy() * q1.threeProduct( phi1 );
const double aux2_2 = itermLL * ( mf2_ + 4.*z2*z2*t2abs ) * phi2.energy2()
+itermTT * ( ( z2p*z2p + z2m*z2m )*( dot2*dot2 + cross2*cross2 ) )
+itermtt * ( cross2*cross2 - dot2*dot2 )
-itermLT * 4.*z2*( z2p-z2m ) * phi2.energy() * q2.threeProduct( phi2 );
//=================================================================
// convention of matrix element as in our kt-factorization
// for heavy flavours
//=================================================================
const double amat2_1 = aux2_1*2.*z1*q1.pt2()/( q1.pt2()*q2.pt2() ),
amat2_2 = aux2_2*2.*z2*q2.pt2()/( q1.pt2()*q2.pt2() );
//=================================================================
// symmetrization
//=================================================================
CG_DEBUG_LOOP( "PPtoFF:offShell" )
<< "aux2(1/2) = " << aux2_1 << " / " << aux2_2 << "\n\t"
<< "amat2(1/2), amat2 = " << amat2_1 << " / " << amat2_2 << " / " << ( 0.5*( imat1*amat2_1 + imat2*amat2_2 ) ) << ".";
return 0.5*( imat1*amat2_1 + imat2*amat2_2 );
}
// register process and define aliases
REGISTER_PROCESS( pptoll, PPtoFF )
REGISTER_PROCESS( pptoff, PPtoFF )
}
}
-
diff --git a/CepGen/Processes/PPtoFF.h b/CepGen/Processes/PPtoFF.h
index 52d32d1..d5e696c 100644
--- a/CepGen/Processes/PPtoFF.h
+++ b/CepGen/Processes/PPtoFF.h
@@ -1,55 +1,54 @@
#ifndef CepGen_Processes_PPtoFF_h
#define CepGen_Processes_PPtoFF_h
#include "CepGen/Processes/GenericKTProcess.h"
#include "CepGen/Core/ParametersList.h"
namespace CepGen
{
- namespace Process
+ namespace process
{
/// Compute the matrix element for a CE \f$\gamma\gamma\rightarrow f\bar f\f$ process using \f$k_T\f$-factorization approach
class PPtoFF : public GenericKTProcess
{
public:
PPtoFF( const ParametersList& params = ParametersList() );
ProcessPtr clone( const ParametersList& params ) const override { return ProcessPtr( new PPtoFF( params ) ); }
private:
void preparePhaseSpace() override;
/// \note IncQQbar in pptoll
double computeKTFactorisedMatrixElement() override;
void fillCentralParticlesKinematics() override;
/// Rapidity range for the outgoing fermions
double onShellME( double shat, double that, double uhat ) const;
double offShellME( double, double, double, double, double, double, const Particle::Momentum&, const Particle::Momentum& ) const;
/// PDG id of the fermion pair produced
int pair_;
int method_;
Limits rap_limits_;
/// Rapidity of the first outgoing fermion
double y1_;
/// Rapidity of the first outgoing fermion
double y2_;
/// Transverse momentum difference for the two outgoing fermions
double pt_diff_;
/// Azimuthal angle difference for the two outgoing fermions
double phi_pt_diff_;
double mf_, mf2_, qf_;
unsigned short colf_;
/// First outgoing fermion's momentum
Particle::Momentum p_f1_;
/// Second outgoing fermion's momentum
Particle::Momentum p_f2_;
};
}
}
#endif
-
diff --git a/CepGen/Processes/PPtoWW.cpp b/CepGen/Processes/PPtoWW.cpp
index f2f5871..4ccf458 100644
--- a/CepGen/Processes/PPtoWW.cpp
+++ b/CepGen/Processes/PPtoWW.cpp
@@ -1,394 +1,393 @@
#include "CepGen/Processes/PPtoWW.h"
#include "CepGen/Event/Event.h"
#include "CepGen/Physics/Constants.h"
#include "CepGen/Physics/FormFactors.h"
#include "CepGen/Physics/PDG.h"
#include "CepGen/Core/Exception.h"
#include <assert.h>
#include "CepGen/Processes/ProcessesHandler.h"
namespace CepGen
{
- namespace Process
+ namespace process
{
const double PPtoWW::mw_ = ParticleProperties::mass( PDG::W );
const double PPtoWW::mw2_ = PPtoWW::mw_*PPtoWW::mw_;
PPtoWW::PPtoWW( const ParametersList& params ) :
GenericKTProcess( params, "pptoww", "ɣɣ → W⁺W¯", { { PDG::photon, PDG::photon } }, { PDG::W, PDG::W } ),
method_( params.get<int>( "method", 1 ) ),
pol_state_( (Polarisation)params.get<int>( "polarisationStates", 0 ) ),
y1_( 0. ), y2_( 0. ), pt_diff_( 0. ), phi_pt_diff_( 0. )
{}
void
PPtoWW::preparePhaseSpace()
{
registerVariable( y1_, Mapping::linear, cuts_.cuts.central.rapidity_single, { -6., 6. }, "First outgoing W rapidity" );
registerVariable( y2_, Mapping::linear, cuts_.cuts.central.rapidity_single, { -6., 6. }, "Second outgoing W rapidity" );
registerVariable( pt_diff_, Mapping::linear, cuts_.cuts.central.pt_diff, { 0., 500. }, "Ws transverse momentum difference" );
registerVariable( phi_pt_diff_, Mapping::linear, cuts_.cuts.central.phi_pt_diff, { 0., 2.*M_PI }, "Ws azimuthal angle difference" );
switch ( pol_state_ ) {
case Polarisation::LL: pol_w1_ = pol_w2_ = { 0 }; break;
case Polarisation::LT: pol_w1_ = { 0 }; pol_w2_ = { -1, 1 }; break;
case Polarisation::TL: pol_w1_ = { -1, 1 }; pol_w2_ = { 0 }; break;
case Polarisation::TT: pol_w1_ = pol_w2_ = { -1, 1 }; break;
default:
case Polarisation::full: pol_w1_ = pol_w2_ = { -1, 0, 1 }; break;
}
CG_DEBUG( "PPtoWW:mode" )
<< "matrix element computation method: " << method_ << ".";
}
double
PPtoWW::computeKTFactorisedMatrixElement()
{
//=================================================================
// 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_ );
CG_DEBUG_LOOP( "PPtoWW:qt" )
<< "q1t(x/y) = " << q1tx << " / " << q1ty << "\n\t"
<< "q2t(x/y) = " << 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 = ( ptsumx+ptdiffx )*0.5, pt1y = ( ptsumy+ptdiffy )*0.5, pt1 = std::hypot( pt1x, pt1y ),
pt2x = ( ptsumx-ptdiffx )*0.5, pt2y = ( ptsumy-ptdiffy )*0.5, pt2 = std::hypot( pt2x, pt2y );
if ( cuts_.cuts.central_particles.count( PDG::W ) > 0
&& cuts_.cuts.central_particles.at( PDG::W ).pt_single.valid() ) {
const Limits pt_limits = cuts_.cuts.central_particles.at( PDG::W ).pt_single;
if ( !pt_limits.passes( pt1 ) || !pt_limits.passes( pt2 ) )
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 );
if ( !cuts_.cuts.central.mass_sum.passes( invm ) )
return 0.;
//=================================================================
// a window in transverse momentum difference
//=================================================================
if ( !cuts_.cuts.central.pt_diff.passes( fabs( pt1-pt2 ) ) )
return 0.;
//=================================================================
// a window in rapidity distance
//=================================================================
if ( !cuts_.cuts.central.rapidity_diff.passes( fabs( y1_-y2_ ) ) )
return 0.;
//=================================================================
// auxiliary quantities
//=================================================================
const double alpha1 = amt1/sqs_*exp( y1_ ), beta1 = amt1/sqs_*exp( -y1_ ),
alpha2 = amt2/sqs_*exp( y2_ ), beta2 = amt2/sqs_*exp( -y2_ );
CG_DEBUG_LOOP( "PPtoWW:sudakov" )
<< "Sudakov parameters:\n\t"
<< " alpha1/2 = " << alpha1 << " / " << alpha2 << "\n\t"
<< " beta1/2 = " << beta1 << " / " << beta2 << ".";
const double q1t2 = q1tx*q1tx+q1ty*q1ty, q2t2 = q2tx*q2tx+q2ty*q2ty;
const double x1 = alpha1+alpha2, x2 = beta1+beta2;
const double z1p = alpha1/x1, z1m = alpha2/x1,
z2p = beta1 /x2, z2m = beta2 /x2;
CG_DEBUG_LOOP( "PPtoWW:zeta" )
<< "z(1/2)p = " << z1p << " / " << z2p << "\n\t"
<< "z(1/2)m = " << 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();
CG_DEBUG_LOOP( "PPtoWW:incoming" )
<< "incoming particles: p1: " << ak1z << " / " << ak10 << "\n\t"
<< " p2: " << ak2z << " / " << ak20 << ".";
//=================================================================
// additional conditions for energy-momentum conservation
//=================================================================
const double s1_eff = x1*s_-qt1_*qt1_, s2_eff = x2*s_-qt2_*qt2_;
CG_DEBUG_LOOP( "PPtoWW:central" )
<< "s(1/2)_eff = " << s1_eff << " / " << s2_eff << " GeV^2\n\t"
<< "diboson invariant mass = " << invm << " GeV";
if ( ( cuts_.mode == KinematicsMode::ElasticInelastic
|| cuts_.mode == KinematicsMode::InelasticInelastic )
&& ( sqrt( s1_eff ) <= ( MY_+invm ) ) )
return 0.;
if ( ( cuts_.mode == KinematicsMode::InelasticElastic
|| cuts_.mode == KinematicsMode::InelasticInelastic )
&& ( sqrt( s2_eff ) <= ( MX_+invm ) ) )
return 0.;
//const double qcaptx = pcaptx, qcapty = pcapty;
//=================================================================
// four-momenta of the outgoing protons (or remnants)
//=================================================================
const double px_plus = ( 1.-x1 )*fabs( ak1z )*M_SQRT2,
px_minus = ( MX_*MX_ + q1t2 )*0.5/px_plus;
const double py_minus = ( 1.-x2 )*fabs( ak2z )*M_SQRT2, // warning! sign of pz??
py_plus = ( MY_*MY_ + q2t2 )*0.5/py_minus;
CG_DEBUG_LOOP( "PPtoWW:pxy" )
<< "px± = " << px_plus << " / " << px_minus << "\n\t"
<< "py± = " << py_plus << " / " << py_minus << ".";
PX_ = Particle::Momentum( -q1tx, -q1ty, ( px_plus-px_minus )*M_SQRT1_2, ( px_plus+px_minus )*M_SQRT1_2 );
PY_ = Particle::Momentum( -q2tx, -q2ty, ( py_plus-py_minus )*M_SQRT1_2, ( py_plus+py_minus )*M_SQRT1_2 );
CG_DEBUG_LOOP( "PPtoWW:remnants" )
<< "First remnant: " << PX_ << ", mass = " << PX_.mass() << "\n\t"
<< "Second remnant: " << PY_ << ", mass = " << PY_.mass() << ".";
/*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 /////////////////////
const 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_ ) );
const Particle::Momentum 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_ ) );
//////////////////////////////////////////
CG_DEBUG_LOOP( "PPtoWW:partons" )
<< "First photon*: " << q1 << ", mass2 = " << q1.mass2() << "\n\t"
<< "Second photon*: " << q2 << ", mass2 = " << 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 );
CG_DEBUG_LOOP( "PPtoWW:central" )
<< "First W: " << p_w1_ << ", mass = " << p_w1_.mass() << "\n\t"
<< "Second W: " << p_w2_ << ", mass = " << 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();
const double uhat1 = ( q1-p_w2_ ).mass2(), uhat2 = ( q2-p_w1_ ).mass2();
CG_DEBUG_LOOP( "PPtoWW" )
<< "that(1/2) = " << that1 << " / " << that2 << "\n\t"
<< "uhat(1/2) = " << uhat1 << " / " << uhat2 << ".";
//const double mll = sqrt( shat );
const double that = 0.5*( that1+that2 ), uhat = 0.5*( uhat1+uhat2 );
//=================================================================
// matrix elements
//=================================================================
double amat2 = 0.;
//=================================================================
// How matrix element is calculated
//=================================================================
CG_DEBUG_LOOP( "PPtoWW" )
<< "matrix element mode: " << method_ << ".";
//=================================================================
// matrix element for gamma gamma --> W^+ W^-
// (Denner+Dittmaier+Schuster)
// (work in collaboration with C. Royon)
//=================================================================
if ( method_ == 0 )
amat2 = onShellME( shat, that, uhat );
//=================================================================
// off-shell Nachtmann formulae
//=================================================================
else if ( method_ == 1 )
amat2 = offShellME( shat, that, uhat, phi_qt1_+phi_qt2_, phi_qt1_-phi_qt2_ );
if ( amat2 <= 0. )
return 0.;
//============================================
// unintegrated photon distributions
//============================================
const std::pair<double,double> fluxes
= GenericKTProcess::incomingFluxes( x1, q1t2, x2, q2t2 );
CG_DEBUG_LOOP( "PPtoWW:fluxes" )
<< "Incoming photon fluxes for (x/kt2) = "
<< "(" << x1 << "/" << q1t2 << "), "
<< "(" << x2 << "/" << q2t2 << "):\n\t"
<< fluxes.first << ", " << fluxes.second << ".";
//=================================================================
// 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 / ( 16.*M_PI*M_PI*( x1*x2*s_ )*( x1*x2*s_ ) )
* fluxes.first*M_1_PI * fluxes.second*M_1_PI * 0.25
* Constants::GeV2toBarn;
/*const double aintegral = amat2 / ( 16.*M_PI*M_PI*x1*x1*x2*x2*s_*s_ )
* fluxes.first*M_1_PI * fluxes.second*M_1_PI
* Constants::GeV2toBarn * 0.25;*/
//=================================================================
return aintegral*qt1_*qt2_*pt_diff_;
//=================================================================
}
void
PPtoWW::fillCentralParticlesKinematics()
{
// randomise the charge of the outgoing leptons
short sign = ( drand() > 0.5 ) ? +1 : -1;
//=================================================================
// first outgoing lepton
//=================================================================
Particle& ow1 = event_->getByRole( Particle::CentralSystem )[0];
ow1.setPdgId( ow1.pdgId(), sign );
ow1.setStatus( Particle::Status::Undecayed );
ow1.setMomentum( p_w1_ );
//=================================================================
// second outgoing lepton
//=================================================================
Particle& ow2 = event_->getByRole( Particle::CentralSystem )[1];
ow2.setPdgId( ow2.pdgId(), -sign );
ow2.setStatus( Particle::Status::Undecayed );
ow2.setMomentum( p_w2_ );
}
double
PPtoWW::onShellME( double shat, double that, double uhat )
{
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.*pow( mw2_-that, 2 )*pow( mw2_-uhat, 2 ) );
const double auxil_gamgam = 1.-term1+term2;
const double beta = sqrt( 1.-4.*mw2_/shat );
return 3.*Constants::alphaEM*Constants::alphaEM*beta / ( 2.*shat ) * auxil_gamgam / ( beta/( 64.*M_PI*M_PI*shat ) );
}
double
PPtoWW::offShellME( double shat, double that, double uhat, double phi_sum, double phi_diff )
{
const double e2 = 4.*M_PI*Constants::alphaEM;
double amat2_0 = 0., amat2_1 = 0., amat2_interf = 0.;
for ( const auto lam3 : pol_w1_ )
for ( const auto lam4 : pol_w2_ ) {
double ampli_pp = amplitudeWW( shat, that, uhat, +1, +1, lam3, lam4 );
double ampli_mm = amplitudeWW( shat, that, uhat, -1, -1, lam3, lam4 );
double ampli_pm = amplitudeWW( shat, that, uhat, +1, -1, lam3, lam4 );
double ampli_mp = amplitudeWW( shat, that, uhat, -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 ) );
}
return e2*e2*( amat2_0+amat2_1+amat2_interf );
}
double
PPtoWW::amplitudeWW( double shat, double that, double uhat, short lam1, short lam2, short lam3, short lam4 )
{
//--- first compute some kinematic variables
const double cos_theta = ( that-uhat ) / shat / sqrt( 1.+1.e-10-4.*mw2_/shat ),
cos_theta2 = cos_theta*cos_theta;
const double sin_theta2 = 1.-cos_theta2,
sin_theta = sqrt( sin_theta2 );
const double beta = sqrt( 1.-4.*mw2_/shat ), beta2 = beta*beta;
const double inv_gamma = sqrt( 1.-beta2 ), gamma = 1./inv_gamma,
gamma2 = gamma*gamma, inv_gamma2 = inv_gamma*inv_gamma;
const double invA = 1./( 1.-beta2*cos_theta2 );
//--- per-helicity amplitude
// longitudinal-longitudinal
if ( lam3 == 0 && lam4 == 0 )
return invA*inv_gamma2*( ( gamma2+1. )*( 1.-lam1*lam2 )*sin_theta2 - ( 1.+lam1*lam2 ) );
// transverse-longitudinal
if ( lam4 == 0 )
return invA*( -M_SQRT2*inv_gamma*( lam1-lam2 )*( 1.+lam1*lam3*cos_theta )*sin_theta );
// longitudinal-transverse
if ( lam3 == 0 )
return invA*( -M_SQRT2*inv_gamma*( lam2-lam1 )*( 1.+lam2*lam4*cos_theta )*sin_theta );
// transverse-transverse
if ( lam3 != 0 && lam4 != 0 )
return -0.5*invA*( 2.*beta*( lam1+lam2 )*( lam3+lam4 )
-inv_gamma2*( 1.+lam3*lam4 )*( 2.*lam1*lam2+( 1.-lam1*lam2 ) * cos_theta2 )
+( 1.+lam1*lam2*lam3*lam4 )*( 3.+lam1*lam2 )
+2.*( lam1-lam2 )*( lam3-lam4 )*cos_theta
+( 1.-lam1*lam2 )*( 1.-lam3*lam4 )*cos_theta2 );
return 0.;
}
// register process and define aliases
REGISTER_PROCESS( pptoww, PPtoWW )
}
}
-
diff --git a/CepGen/Processes/PPtoWW.h b/CepGen/Processes/PPtoWW.h
index 0fb2e68..49a10c1 100644
--- a/CepGen/Processes/PPtoWW.h
+++ b/CepGen/Processes/PPtoWW.h
@@ -1,59 +1,58 @@
#ifndef CepGen_Processes_PPtoWW_h
#define CepGen_Processes_PPtoWW_h
#include "CepGen/Processes/GenericKTProcess.h"
#include "CepGen/Core/ParametersList.h"
namespace CepGen
{
- namespace Process
+ namespace process
{
/// \brief Compute the matrix element for a CE \f$\gamma\gamma\rightarrow W^+W^-\f$ process using \f$k_T\f$-factorization approach
/// \note The full theoretical description of this process definition may be found in \cite Luszczak:2018ntp.
class PPtoWW : public GenericKTProcess
{
public:
PPtoWW( const ParametersList& params = ParametersList() );
ProcessPtr clone( const ParametersList& params ) const override { return ProcessPtr( new PPtoWW( params ) ); }
enum class Polarisation { full = 0, LL = 1, LT = 2, TL = 3, TT = 4 };
private:
static const double mw_, mw2_;
void preparePhaseSpace() override;
double computeKTFactorisedMatrixElement() override;
void fillCentralParticlesKinematics() override;
double amplitudeWW( double shat, double that, double uhat, short lam1, short lam2, short lam3, short lam4 );
double onShellME( double shat, double that, double uhat );
double offShellME( double shat, double that, double uhat, double phi_sum, double phi_diff );
int method_;
Polarisation pol_state_;
std::vector<short> pol_w1_, pol_w2_;
/// Rapidity range for the outgoing W bosons
Limits rap_limits_;
/// Rapidity of the first outgoing W boson
double y1_;
/// Rapidity of the first outgoing W boson
double y2_;
Limits ptdiff_limits_;
/// Transverse momentum difference for the two outgoing W bosons
double pt_diff_;
Limits phi_pt_diff_limits_;
/// Azimuthal angle difference for the two outgoing W bosons
double phi_pt_diff_;
/// First outgoing W boson's momentum
Particle::Momentum p_w1_;
/// Second outgoing W boson's momentum
Particle::Momentum p_w2_;
};
}
}
#endif
-
diff --git a/CepGen/Processes/ProcessesHandler.cpp b/CepGen/Processes/ProcessesHandler.cpp
index b5c8ab4..312f2b8 100644
--- a/CepGen/Processes/ProcessesHandler.cpp
+++ b/CepGen/Processes/ProcessesHandler.cpp
@@ -1,44 +1,44 @@
#include "CepGen/Core/Exception.h"
#include "CepGen/Core/utils.h"
#include "CepGen/Core/ParametersList.h"
#include "CepGen/Processes/PPtoFF.h"
#include "CepGen/Processes/ProcessesHandler.h"
#include <sstream>
namespace CepGen
{
ProcessesHandler&
ProcessesHandler::get()
{
static ProcessesHandler instance;
return instance;
}
void
- ProcessesHandler::registerProcess( const std::string& name, const Process::GenericProcess* proc )
+ ProcessesHandler::registerProcess( const std::string& name, const process::GenericProcess* proc )
{
map_[name].reset( proc );
CG_DEBUG( "ProcessesHandler" ) << "Process name \"" << name << "\" registered in database.";
}
void
ProcessesHandler::dump() const
{
std::ostringstream oss;
for ( const auto& p : map_ )
oss << " '" << p.first << "'";
CG_INFO( "ProcessesHandler:dump" )
<< "List of process(es) handled in the database:" << oss.str();
}
ProcessPtr
ProcessesHandler::build( const std::string& name, const ParametersList& params ) const
{
if ( map_.count( name ) == 0 )
throw CG_FATAL( "ProcessesHandler:build" )
<< "Failed to retrieve a process with name \"" << name << "\"!";
return map_.at( name )->clone( params );
}
}
diff --git a/CepGen/Processes/ProcessesHandler.h b/CepGen/Processes/ProcessesHandler.h
index bef448b..f7c23b7 100644
--- a/CepGen/Processes/ProcessesHandler.h
+++ b/CepGen/Processes/ProcessesHandler.h
@@ -1,40 +1,39 @@
#ifndef CepGen_Processes_ProcessesHandler_h
#define CepGen_Processes_ProcessesHandler_h
#include "CepGen/Processes/GenericProcess.h"
#include <unordered_map>
#include <memory>
#define BUILDERNM( obj ) obj ## Builder
#define STRINGIFY( name ) #name
#define REGISTER_PROCESS( name, obj ) \
struct BUILDERNM( name ) { \
BUILDERNM( name )() { CepGen::ProcessesHandler::get().registerProcess( STRINGIFY( name ), new obj ); } }; \
static BUILDERNM( name ) g ## name;
namespace CepGen
{
class ParametersList;
- //namespace Process { class GenericProcess; }
+ //namespace process { class GenericProcess; }
class ProcessesHandler
{
public:
static ProcessesHandler& get();
~ProcessesHandler() = default;
- void registerProcess( const std::string& name, const CepGen::Process::GenericProcess* );
+ void registerProcess( const std::string& name, const CepGen::process::GenericProcess* );
ProcessPtr build( const std::string& name, const ParametersList& ) const;
void dump() const;
private:
explicit ProcessesHandler() = default;
- std::unordered_map<std::string, std::unique_ptr<const Process::GenericProcess> > map_;
+ std::unordered_map<std::string, std::unique_ptr<const process::GenericProcess> > map_;
public:
ProcessesHandler( const ProcessesHandler& ) = delete;
void operator=( const ProcessesHandler& ) = delete;
};
}
#endif
-
diff --git a/test/cepgen.cpp b/test/cepgen.cpp
index f036e0a..9cd8999 100644
--- a/test/cepgen.cpp
+++ b/test/cepgen.cpp
@@ -1,74 +1,74 @@
//--- steering cards
#include "CepGen/Cards/PythonHandler.h"
#include "CepGen/Cards/LpairHandler.h"
#include "CepGen/Generator.h"
#include "CepGen/Core/Exception.h"
//--- necessary include to build the default run
#include "CepGen/Processes/GamGamLL.h"
#include "CepGen/Physics/PDG.h"
#include "abort.h"
#include <iostream>
using namespace std;
/**
* Main caller for this MC generator.
* * loads the configuration files' variables if passed as an argument,
* or a default LPAIR-like configuration,
* * launches the cross-section computation and the events generation.
* \author Laurent Forthomme <laurent.forthomme@cern.ch>
*/
int main( int argc, char* argv[] ) {
//--- first start by defining the generator object
CepGen::Generator gen;
if ( argc < 2 ) {
CG_INFO( "main" ) << "No config file provided. Setting the default parameters.";
//--- default run: LPAIR elastic ɣɣ → µ⁺µ¯ at 13 TeV
CepGen::ParametersList pgen;
pgen.set<int>( "pair", (int)CepGen::PDG::muon );
- gen.parameters->setProcess( new CepGen::Process::GamGamLL( pgen ) );
+ gen.parameters->setProcess( new CepGen::process::GamGamLL( pgen ) );
gen.parameters->kinematics.mode = CepGen::KinematicsMode::ElasticElastic;
gen.parameters->kinematics.cuts.central.pt_single.min() = 15.;
gen.parameters->kinematics.cuts.central.eta_single = { -2.5, 2.5 };
gen.parameters->generation.enabled = true;
gen.parameters->generation.maxgen = 1e3;
}
else {
CG_INFO( "main" ) << "Reading config file stored in " << argv[1] << ".";
const std::string extension = CepGen::Cards::Handler::getExtension( argv[1] );
if ( extension == "card" )
gen.setParameters( CepGen::Cards::LpairHandler( argv[1] ).parameters() );
#ifdef PYTHON
else if ( extension == "py" )
gen.setParameters( CepGen::Cards::PythonHandler( argv[1] ).parameters() );
#endif
else
throw CG_FATAL( "main" ) << "Unrecognized steering card extension: ." << extension << "!";
}
//--- list all parameters
CG_INFO( "main" ) << gen.parameters.get();
AbortHandler ctrl_c;
try {
//--- let there be a cross-section...
double xsec = 0., err = 0.;
gen.computeXsection( xsec, err );
if ( gen.parameters->generation.enabled )
//--- events generation starts here
// (one may use a callback function)
gen.generate();
} catch ( const CepGen::RunAbortedException& e ) {
CG_INFO( "main" ) << "Run aborted!";
}
return 0;
}