diff --git a/CepGen/Core/Exception.h b/CepGen/Core/Exception.h
index 800b4ab..eb610ae 100644
--- a/CepGen/Core/Exception.h
+++ b/CepGen/Core/Exception.h
@@ -1,202 +1,202 @@
#ifndef CepGen_Core_Exception_h
#define CepGen_Core_Exception_h
#include <csignal>
#include "CepGen/Core/Logger.h"
namespace cepgen
{
/// A generic exception type
/// \author Laurent Forthomme <laurent.forthomme@cern.ch>
/// \date 27 Mar 2015
struct Exception
{
explicit inline Exception() = default;
virtual ~Exception() noexcept = default;
/// Enumeration of exception severities
enum class Type {
undefined = -1, ///< Irregular exception
debug, ///< Debugging information to be enabled
verbatim, ///< Raw information
info, ///< Prettified information
warning, ///< Casual non-stopping warning
error, ///< General non-stopping error
fatal ///< Critical and stopping error
};
/// Stream operator (null and void)
template<class T> Exception& operator<<( const T& ) { return *this; }
/// Dump the full exception information in a given output stream
/// \param[inout] os the output stream where the information is dumped
virtual void dump( std::ostream& os = *utils::Logger::get().output ) const = 0;
/// Exception message
virtual std::string message() const = 0;
};
/// A simple exception handler
/// \date 24 Mar 2015
class LoggedException : public Exception
{
public:
/// Generic constructor
/// \param[in] module exception classifier
/// \param[in] type exception type
/// \param[in] id exception code (useful for logging)
explicit inline LoggedException( const char* module = "", Type type = Type::undefined, short id = 0 ) :
module_( module ), type_( type ), error_num_( id ) {}
/// Generic constructor
/// \param[in] from method invoking the exception
/// \param[in] module exception classifier
/// \param[in] type exception type
/// \param[in] id exception code (useful for logging)
explicit inline LoggedException( const char* from, const char* module, Type type = Type::undefined, short id = 0 ) :
from_( from ), module_( module ), type_( type ), error_num_( id ) {}
/// Copy constructor
inline LoggedException( const LoggedException& rhs ) :
from_( rhs.from_ ), module_( rhs.module_ ), message_( rhs.message_.str() ),
type_( rhs.type_ ), error_num_( rhs.error_num_ ) {}
/// Default destructor (potentially killing the process)
inline ~LoggedException() noexcept override {
if ( type_ != Type::undefined )
dump();
// we stop this process' execution on fatal exception
if ( type_ == Type::fatal && raise( SIGINT ) != 0 )
exit( 0 );
}
//----- Overloaded stream operators
/// Generic templated message feeder operator
template<typename T>
inline friend const LoggedException& operator<<( const LoggedException& exc, T var ) {
LoggedException& nc_except = const_cast<LoggedException&>( exc );
nc_except.message_ << var;
return exc;
}
/// Pipe modifier operator
inline friend const LoggedException& operator<<( const LoggedException& exc, std::ios_base&( *f )( std::ios_base& ) ) {
LoggedException& nc_except = const_cast<LoggedException&>( exc );
f( nc_except.message_ );
return exc;
}
inline std::string message() const override {
return message_.str();
}
/// Extract the origin of the exception
inline std::string from() const { return from_; }
/// Extract the exception code
inline int errorNumber() const { return error_num_; }
/// Extract the exception type
inline Type type() const { return type_; }
/// Extract a human-readable (and colourified) version of the exception type
inline std::string typeString() const {
switch ( type() ) {
case Type::warning: return "\033[34;1mWarning\033[0m";
case Type::info: return "\033[32;1mInfo.\033[0m";
case Type::debug: return "\033[33;1mDebug\033[0m";
case Type::error: return "\033[31;1mError\033[0m";
case Type::fatal: return "\033[31;1mFatal\033[0m";
case Type::undefined: default: return "\33[7;1mUndefined\033[0m";
}
}
inline void dump( std::ostream& os = *utils::Logger::get().output ) const override {
os << fullMessage() << std::endl;
}
/// Extract a one-line summary of the exception
inline std::string shortMessage() const {
std::ostringstream os;
os << "[" << typeString() << "]";
- if ( type_ == Type::warning )
+ if ( type_ == Type::warning || type_ == Type::debug )
os << " \033[30;4m" << from_ << "\033[0m\n";
os << "\t" << message_.str();
return os.str();
}
private:
inline static char* now() {
static char buffer[10];
time_t rawtime;
time( &rawtime );
struct tm* timeinfo = localtime( &rawtime );
strftime( buffer, 10, "%H:%M:%S", timeinfo );
return buffer;
}
/// Extract a full exception message
inline std::string fullMessage() const {
switch ( type_ ) {
case Type::info:
case Type::debug:
case Type::warning:
return shortMessage();
case Type::verbatim:
return message_.str();
default: {
std::ostringstream os;
os << "================================== Exception ========================= " << now()
<< "\n Class: " << typeString() << std::endl;
if ( !from_.empty() )
os << " Raised by: " << from_ << "\n"
<< " " << message_.str() << std::endl;
if ( errorNumber() != 0 )
os << "-------------------------------------------------------------------------------"
<< "\n Error #" << error_num_ << std::endl;
os << "===============================================================================";
return os.str();
}
}
}
/// Origin of the exception
std::string from_;
/// Exception classificator
std::string module_;
/// Message to throw
std::ostringstream message_;
/// Exception type
Type type_;
/// Integer exception number
short error_num_;
};
/// Placeholder for debugging messages if logging threshold is not reached
/// \date Apr 2018
struct NullStream : Exception
{
using Exception::Exception;
/// Empty constructor
inline NullStream() {}
/// Empty constructor
inline NullStream( const LoggedException& ) {}
void dump( std::ostream& os = *utils::Logger::get().output ) const override {}
std::string message() const override { return ""; }
};
}
#define CG_LOG( mod ) \
( !CG_LOG_MATCH( mod, information ) ) \
? cepgen::NullStream() \
: cepgen::LoggedException( __PRETTY_FUNCTION__, mod, cepgen::Exception::Type::verbatim )
#define CG_INFO( mod ) \
( !CG_LOG_MATCH( mod, information ) ) \
? cepgen::NullStream() \
: cepgen::LoggedException( __PRETTY_FUNCTION__, mod, cepgen::Exception::Type::info )
#define CG_DEBUG( mod ) \
( !CG_LOG_MATCH( mod, debug ) ) \
? cepgen::NullStream() \
: cepgen::LoggedException( __PRETTY_FUNCTION__, mod, cepgen::Exception::Type::debug )
#define CG_DEBUG_LOOP( mod ) \
( !CG_LOG_MATCH( mod, debugInsideLoop ) ) \
? cepgen::NullStream() \
: cepgen::LoggedException( __PRETTY_FUNCTION__, mod, cepgen::Exception::Type::debug )
#define CG_WARNING( mod ) \
( !CG_LOG_MATCH( mod, warning ) ) \
? cepgen::NullStream() \
: cepgen::LoggedException( __PRETTY_FUNCTION__, mod, cepgen::Exception::Type::warning )
#define CG_ERROR( mod ) \
cepgen::LoggedException( __PRETTY_FUNCTION__, mod, cepgen::Exception::Type::error )
#define CG_FATAL( mod ) \
cepgen::LoggedException( __PRETTY_FUNCTION__, mod, cepgen::Exception::Type::fatal )
#endif
diff --git a/CepGen/Core/utils.h b/CepGen/Core/utils.h
index bb63dee..e94fac1 100644
--- a/CepGen/Core/utils.h
+++ b/CepGen/Core/utils.h
@@ -1,50 +1,50 @@
#ifndef CepGen_Core_utils_h
#define CepGen_Core_utils_h
#include <string>
#include <vector>
#include <numeric>
namespace cepgen
{
/// Format a string using a printf style format descriptor.
std::string Form( const std::string fmt, ... );
/// Human-readable boolean printout
inline const char* yesno( const bool& test ) { return ( test ) ? "\033[32;1myes\033[0m" : "\033[31;1mno\033[0m"; }
//inline const char* boldify( const char* str ) { const std::string out = std::string( "\033[33;1m" ) + std::string( str ) + std::string( "\033[0m" ); return out.c_str(); }
/// Boldify a string for TTY-type output streams
inline std::string boldify( const std::string& str ) { return Form( "\033[1m%s\033[0m", str.c_str() ); }
/// Boldify a string for TTY-type output streams
inline std::string boldify( const char* str ) { return boldify( std::string( str ) ); }
/// Boldify a double floating point number for TTY-type output streams
inline std::string boldify( const double& dbl ) { return boldify( Form("%.2f", dbl ) ); }
/// Boldify an integer for TTY-type output streams
inline std::string boldify( const int& i ) { return boldify( Form("% d", i ) ); }
/// Boldify an unsigned integer for TTY-type output streams
inline std::string boldify( const unsigned int& ui ) { return boldify( Form("%d", ui ) ); }
/// Boldify an unsigned long integer for TTY-type output streams
inline std::string boldify( const unsigned long& ui ) { return boldify( Form("%lu", ui ) ); }
/// TTY-type enumeration of colours
enum class Colour { gray = 30, red = 31, green = 32, yellow = 33, blue = 34, purple = 35 };
/// Colourise a string for TTY-type output streams
inline std::string colourise( const std::string& str, const Colour& col ) { return Form( "\033[%d%s\033[0m", (int)col, str.c_str() ); }
/// Replace all occurences of a text by another
size_t replace_all( std::string& str, const std::string& from, const std::string& to );
namespace utils
{
/// Add a trailing "s" when needed
inline const char* s( unsigned short num ) { return ( num > 1 ) ? "s" : ""; }
/// Helper to print a vector
template<class T> std::string repr( const std::vector<T>& vec, const std::string& sep = "," ) {
return std::accumulate( std::next( vec.begin() ), vec.end(),
std::to_string( *vec.begin() ), [&sep]( std::string str, T xv ) {
return std::move( str )+sep+std::to_string( xv );
} );
}
}
}
/// Provide a random number generated along a uniform distribution between 0 and 1
-#define drand() static_cast<double>( rand()/RAND_MAX )
+#define drand() (double)rand()/RAND_MAX
#endif
diff --git a/CepGen/Event/Momentum.cpp b/CepGen/Event/Momentum.cpp
index f94fc60..a9ca937 100644
--- a/CepGen/Event/Momentum.cpp
+++ b/CepGen/Event/Momentum.cpp
@@ -1,387 +1,392 @@
#include "CepGen/Event/Particle.h"
#include "CepGen/Core/Exception.h"
#include "CepGen/Core/utils.h"
#include <math.h>
#include <iomanip>
namespace cepgen
{
using Momentum = Particle::Momentum;
//----- Particle momentum methods
Momentum::Momentum() :
px_( 0. ), py_( 0. ), pz_( 0. ), p_( 0. ), energy_( 0. )
{}
Momentum::Momentum( double x, double y, double z, double t ) :
px_( x ), py_( y ), pz_( z ), energy_( t )
{
computeP();
}
Momentum::Momentum( double* p ) :
px_( p[0] ), py_( p[1] ), pz_( p[2] ), energy_( p[3] )
{
computeP();
}
//--- static constructors
Momentum
Momentum::fromPtEtaPhi( double pt, double eta, double phi, double e )
{
const double px = pt*cos( phi ),
py = pt*sin( phi ),
pz = pt*sinh( eta );
return Momentum( px, py, pz, e );
}
Momentum
Momentum::fromPThetaPhi( double p, double theta, double phi, double e )
{
const double px = p*sin( theta )*cos( phi ),
py = p*sin( theta )*sin( phi ),
pz = p*cos( theta );
return Momentum( px, py, pz, e );
}
Momentum
Momentum::fromPxPyPzE( double px, double py, double pz, double e )
{
return Momentum( px, py, pz, e );
}
Momentum
Momentum::fromPxPyPzM( double px, double py, double pz, double m )
{
Momentum mom( px, py, pz );
mom.setMass( m );
return mom;
}
Momentum
Momentum::fromPxPyYM( double px, double py, double rap, double m )
{
const double pt = std::hypot( px, py ), et = std::hypot( pt, m );
return Momentum( px, py, et*sinh( rap ), et*cosh( rap ) );
}
//--- arithmetic operators
Momentum&
Momentum::operator+=( const Momentum& mom )
{
px_ += mom.px_;
py_ += mom.py_;
pz_ += mom.pz_;
energy_ += mom.energy_;
computeP();
return *this;
}
Momentum&
Momentum::operator-=( const Momentum& mom )
{
px_ -= mom.px_;
py_ -= mom.py_;
pz_ -= mom.pz_;
energy_ -= mom.energy_;
computeP();
return *this;
}
bool
Momentum::operator==( const Momentum& mom ) const
{
return ( px_ == mom.px_
&& py_ == mom.py_
&& pz_ == mom.pz_
&& energy_ == mom.energy_ );
}
double
Momentum::threeProduct( const Momentum& mom ) const
{
CG_DEBUG_LOOP( "Momentum" )
<< " (" << px_ << ", " << py_ << ", " << pz_ << ")\n\t"
<< "* (" << mom.px_ << ", " << mom.py_ << ", " << mom.pz_ << ")\n\t"
<< "= " << px_*mom.px_+py_*mom.py_+pz_*mom.pz_;
return px_*mom.px_+py_*mom.py_+pz_*mom.pz_;
}
double
Momentum::fourProduct( const Momentum& mom ) const
{
CG_DEBUG_LOOP( "Momentum" )
<< " (" << px_ << ", " << py_ << ", " << pz_ << ", " << energy_ << ")\n\t"
<< "* (" << mom.px_ << ", " << mom.py_ << ", " << mom.pz_ << ", " << mom.energy_ << ")\n\t"
<< "= " << energy_*mom.energy_-threeProduct(mom);
return energy_*mom.energy_-threeProduct(mom);
}
double
Momentum::crossProduct( const Momentum& mom ) const
{
return px_*mom.py_-py_*mom.px_;
}
double
Momentum::operator*=( const Momentum& mom )
{
return threeProduct( mom );
}
Momentum&
Momentum::operator*=( double c )
{
px_ *= c;
py_ *= c;
pz_ *= c;
energy_ *= c;
computeP();
return *this;
}
Momentum
operator*( const Momentum& mom, double c )
{
Momentum out = mom;
out *= c;
return out;
}
Momentum
operator*( double c, const Momentum& mom )
{
Momentum out = mom;
out *= c;
return out;
}
Momentum
operator+( const Momentum& mom1, const Momentum& mom2 )
{
Momentum out = mom1;
out += mom2;
return out;
}
Momentum
operator-( const Momentum& mom1, const Momentum& mom2 )
{
Momentum out = mom1;
out -= mom2;
return out;
}
Momentum
operator-( const Momentum& mom )
{
return Momentum()-mom;
}
double
operator*( const Momentum& mom1, const Momentum& mom2 )
{
Momentum tmp = mom1;
return tmp.threeProduct( mom2 );
}
//--- various setters
void
Momentum::setMass2( double m2 )
{
energy_ = sqrt( p2()+m2 );
}
void
Momentum::setP( double px, double py, double pz, double e )
{
setP( px, py, pz );
setEnergy( e );
}
void
Momentum::setP( double px, double py, double pz )
{
px_ = px;
py_ = py;
pz_ = pz;
computeP();
}
void
Momentum::computeP()
{
p_ = std::hypot( pt(), pz_ );
}
void
Momentum::truncate( double tolerance )
{
if ( px_ <= tolerance )
px_ = 0.;
if ( py_ <= tolerance )
py_ = 0.;
if ( pz_ <= tolerance )
pz_ = 0.;
if ( energy_ <= tolerance )
energy_ = 0.;
computeP();
}
//--- various getters
double
Momentum::operator[]( unsigned int i ) const
{
switch ( i ) {
case 0: return px_;
case 1: return py_;
case 2: return pz_;
case 3: return energy_;
default:
throw CG_FATAL( "Momentum" ) << "Failed to retrieve the component " << i << "!";
}
}
double&
Momentum::operator[]( const unsigned int i )
{
switch ( i ) {
case 0: return px_;
case 1: return py_;
case 2: return pz_;
case 3: return energy_;
default:
throw CG_FATAL( "Momentum" ) << "Failed to retrieve the component " << i << "!";
}
}
const std::vector<double>
Momentum::pVector() const
{
return std::vector<double>{ px(), py(), pz(), energy(), mass() };
}
double
Momentum::mass() const
{
if ( mass2() >= 0. )
return sqrt( mass2() );
return -sqrt( -mass2() );
}
double
Momentum::theta() const
{
return atan2( pt(), pz() );
}
double
Momentum::phi() const
{
return atan2( py(), px() );
}
double
Momentum::pt() const
{
return std::hypot( px_, py_ );
}
double
Momentum::pt2() const
{
return px_*px_+py_*py_;
}
double
Momentum::eta() const
{
const int sign = ( pz()/fabs( pz() ) );
return ( pt() != 0. )
? log( ( p()+fabs( pz() ) )/pt() )*sign
: 9999.*sign;
}
double
Momentum::rapidity() const
{
const int sign = ( pz()/fabs( pz() ) );
return ( energy() >= 0. )
? log( ( energy()+pz() )/( energy()-pz() ) )*0.5
: 999.*sign;
}
//--- boosts/rotations
Momentum&
Momentum::betaGammaBoost( double gamma, double betagamma )
{
- if ( gamma == 1. && betagamma == 0. ) return *this; // trivial case
+ if ( gamma == 1. && betagamma == 0. )
+ return *this; // trivial case
+
const double pz = pz_, e = energy_;
- pz_ = gamma*pz+betagamma*e;
+ pz_ = gamma*pz+betagamma*e;
energy_ = gamma*e +betagamma*pz;
computeP();
return *this;
}
Momentum&
- Momentum::lorentzBoost( const Momentum& p )
+ Momentum::lorentzBoost( const Momentum& p2 )
{
//--- do not boost on a system at rest
- if ( p.p() == 0. )
+ if ( p2.p() == 0. )
return *this;
- const double m = p.mass();
- const double pf4 = fourProduct( p )/m;
- const double fn = ( pf4+energy_ )/( p.energy_+m );
- *this -= p*fn;
+ const double m = p2.mass();
+ const double pf4 = ( px_*p2.px_+py_*p2.py_+pz_*p2.pz_+energy_*p2.energy_ )/m;
+ const double fn = ( pf4+energy_ )/( p2.energy_+m );
+ *this += p2*fn;
energy_ = pf4;
return *this;
}
Momentum&
Momentum::rotatePhi( double phi, double sign )
{
- const double px = px_*cos( phi )+py_*sin( phi )*sign,
- py =-px_*sin( phi )+py_*cos( phi )*sign;
+ const double sphi = sin( phi ), cphi = cos( phi );
+ const double px = px_*cphi + sign*py_*sphi,
+ py = -px_*sphi + sign*py_*cphi;
px_ = px;
py_ = py;
return *this;
}
Momentum&
Momentum::rotateThetaPhi( double theta, double phi )
{
+ const double ctheta = cos( theta ), stheta = sin( theta );
+ const double cphi = cos( phi ), sphi = sin( phi );
double rotmtx[3][3], mom[3]; //FIXME check this! cos(phi)->-sin(phi) & sin(phi)->cos(phi) --> phi->phi+pi/2 ?
- rotmtx[0][0] = -sin( phi ); rotmtx[0][1] = -cos( theta )*cos( phi ); rotmtx[0][2] = sin( theta )*cos( phi );
- rotmtx[1][0] = cos( phi ); rotmtx[1][1] = -cos( theta )*sin( phi ); rotmtx[1][2] = sin( theta )*sin( phi );
- rotmtx[2][0] = 0.; rotmtx[2][1] = sin( theta ); rotmtx[2][2] = cos( theta );
+ rotmtx[0][0] = -sphi; rotmtx[0][1] = -ctheta*cphi; rotmtx[0][2] = stheta*cphi;
+ rotmtx[1][0] = cphi; rotmtx[1][1] = -ctheta*sphi; rotmtx[1][2] = stheta*sphi;
+ rotmtx[2][0] = 0.; rotmtx[2][1] = stheta; rotmtx[2][2] = ctheta;
for ( unsigned short i = 0; i < 3; ++i ) {
mom[i] = 0.;
for ( unsigned short j = 0; j < 3; ++j )
mom[i] += rotmtx[i][j]*operator[]( j );
}
setP( mom[0], mom[1], mom[2] );
return *this;
}
//--- printout
std::ostream&
operator<<( std::ostream& os, const Momentum& mom )
{
return os << "(" << std::fixed
<< std::setw( 9 ) << mom.energy_ << "|"
<< std::setw( 9 ) << mom.px_ << " "
<< std::setw( 9 ) << mom.py_ << " "
<< std::setw( 9 ) << mom.pz_ << ")" << std::defaultfloat;
}
}
diff --git a/CepGen/Hadronisers/Pythia6Hadroniser.cpp b/CepGen/Hadronisers/Pythia6Hadroniser.cpp
index 4885cca..de6bfb3 100644
--- a/CepGen/Hadronisers/Pythia6Hadroniser.cpp
+++ b/CepGen/Hadronisers/Pythia6Hadroniser.cpp
@@ -1,345 +1,344 @@
#include "CepGen/Hadronisers/GenericHadroniser.h"
#include "CepGen/Hadronisers/HadronisersHandler.h"
#include "CepGen/Event/Event.h"
#include "CepGen/Event/Particle.h"
#include "CepGen/Physics/PDG.h"
#include "CepGen/Core/ParametersList.h" //FIXME
#include "CepGen/Core/Exception.h"
#include "CepGen/Core/utils.h"
#include <algorithm>
/// Maximal number of characters to fetch for the particle's name
#define NAME_CHR 16
extern "C"
{
/// Get the particle's mass in GeV from the Pythia6 module
extern double pymass_( int& );
/// Get the resonant particle's width in GeV from the Pythia6 module
//extern double pywidt_( int& );
/// Launch the Pythia6 fragmentation
extern void pyexec_();
/// Set a parameter value to the Pythia6 module
extern void pygive_( const char*, int );
extern void pyckbd_();
/// List all the particles in the event in a human-readable format
extern void pylist_( int& );
/// Join two coloured particles in a colour singlet
extern void pyjoin_( int&, int& );
/// Get a particle's human-readable name from the Pythia6 module
extern void pyname_( int&, char*, int );
/// Get kinematic information on a particle from the Pythia6 module
extern double pyp_( int&, int& );
/// Store one parton/particle in the PYJETS common block
extern void py1ent_( int&, int&, double&, double&, double& );
void pystop_() {}
/// Particles content of the event
extern struct
{
/// Number of particles in the event
int n;
int npad;
/// Particles' general information (status, PDG id, mother, daughter 1, daughter 2)
int k[5][4000];
/// Particles' kinematics, in GeV (px, py, pz, E, M)
double p[5][4000];
/// Primary vertex for the particles
double v[5][4000];
} pyjets_;
}
namespace cepgen
{
namespace hadr
{
/**
* Full interface to the Pythia6 @cite Sjostrand:2006za 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 Pythia6 hadronisation algorithm
*/
class Pythia6Hadroniser : public GenericHadroniser
{
public:
Pythia6Hadroniser( const ParametersList& );
void setParameters( const Parameters& ) override {}
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 hadronise( const Particle* );
private:
static constexpr unsigned short MAX_PART_STRING = 3;
static constexpr unsigned short MAX_STRING_EVENT = 2;
inline static double pymass(int pdgid_) { return pymass_(pdgid_); }
//inline static double pywidt(int pdgid_) { return pywidt_(pdgid_); }
inline static void pyexec() { pyexec_(); }
inline static void pyckbd() { pyckbd_(); }
inline static void pygive( const std::string& line ) { pygive_( line.c_str(), line.length() ); }
inline static void pylist( int mlist ) { pylist_( mlist ); }
inline static double pyp( int role, int qty ) { return pyp_( role, qty ); }
//inline static void py1ent( int* kf, double* pe, double theta, double phi ) { int one=1; py1ent_( &one, kf, pe, theta, phi ); }
//inline static void py1ent( int* kf, double* pe, double theta, double phi ) { py1ent_( 1, kf, pe, theta, phi ); }
inline static std::string pyname( int pdgid ) {
char out[NAME_CHR];
std::string s;
pyname_( pdgid, out, NAME_CHR );
s = std::string( out, NAME_CHR );
//s.erase(std::find_if(s.rbegin(), s.rend(), std::not1(std::ptr_fun<int, int>(std::isspace))).base(), s.end());
s.erase( remove( s.begin(), s.end(), ' ' ), s.end() );
return s;
}
/**
* \brief Connect entries with colour flow information
* \param[in] njoin_ Number of particles to join in the colour flow
* \param[in] ijoin_ List of particles unique identifier to join in the colour flow
*/
inline static void pyjoin( int njoin, int ijoin[2] ) { return pyjoin_( njoin, *ijoin ); }
bool prepareHadronisation( Event& );
+ std::pair<short,short> pickPartonsContent() const;
struct EventProperties
{
unsigned int str_in_evt = 0;
unsigned int num_part_in_str[MAX_STRING_EVENT] = { 0 };
int jlrole[MAX_STRING_EVENT] = { -1 };
int jlpsf[MAX_STRING_EVENT][MAX_PART_STRING] = { 0 };
};
EventProperties fillParticles( const Event& ) const;
};
Pythia6Hadroniser::Pythia6Hadroniser( const ParametersList& plist ) :
GenericHadroniser( plist, "pythia6" )
{}
void
Pythia6Hadroniser::readString( const char* param )
{
pygive( param );
}
bool
Pythia6Hadroniser::run( Event& ev, double& weight, bool full )
{
weight = 1.;
- if ( full )
+ //if ( full )
prepareHadronisation( ev );
if ( utils::Logger::get().level >= utils::Logger::Level::debug ) {
CG_DEBUG( "Pythia6Hadroniser" ) << "Dump of the event before the hadronisation:";
ev.dump();
}
//--- fill Pythia 6 common blocks
auto prop = fillParticles( ev );
CG_DEBUG( "Pythia6Hadroniser" )
<< "Passed the string construction stage.\n\t"
<< " " << prop.str_in_evt << " string objects were identified and constructed.";
unsigned int oldnpart = pyjets_.n;
for ( unsigned short i = 0; i < prop.str_in_evt; ++i ) {
if ( prop.num_part_in_str[i] < 2 )
continue;
std::ostringstream dbg;
for ( unsigned short j = 0; j < prop.num_part_in_str[i]; ++j )
if ( prop.jlpsf[i][j] != -1 )
dbg << Form( "\n\t * %2d (pdgId=%4d)", prop.jlpsf[i][j], pyjets_.k[1][prop.jlpsf[i][j]-1] );
CG_INFO( "Pythia6Hadroniser" )
<< "Joining " << prop.num_part_in_str[i] << " particle" << utils::s( prop.num_part_in_str[i] )
<< " in a same string (" << i << ")"
<< dbg.str();
pyjoin( prop.num_part_in_str[i], prop.jlpsf[i] );
}
pyexec();
int criteria = oldnpart+1;
for ( unsigned int i = 0; i < MAX_STRING_EVENT; ++i )
criteria += prop.num_part_in_str[i];
if ( pyjets_.k[1][criteria] == 2212
&& pyjets_.k[0][criteria] == 1 ) {
CG_WARNING( "Pythia6Hadroniser" ) << "System is non-inelastic.";
return false;
}
// We filter the first particles already present in the event
for ( unsigned int p = oldnpart; p < (unsigned int)pyjets_.n; ++p ) {
const Particle::Role role = pyjets_.k[2][p] != 0
? ev[pyjets_.k[2][p]-1].role() // child particle inherits its mother's role
: Particle::Role::UnknownRole;
auto& pa = ev.addParticle( role );
pa.setId( p );
pa.setPdgId( (short)pyjets_.k[1][p] );
pa.setStatus( static_cast<Particle::Status>( pyjets_.k[0][p] ) );
pa.setMomentum( Particle::Momentum( pyjets_.p[0][p], pyjets_.p[1][p], pyjets_.p[2][p], pyjets_.p[3][p] ) );
pa.setMass( pyjets_.p[4][p] );
//pa.setCharge( (float)pyp( p+1, 6 ) );
if ( role != Particle::Role::UnknownRole ) {
auto& moth = ev[pyjets_.k[2][p]-1];
moth.setStatus( role == Particle::Role::CentralSystem
? Particle::Status::Resonance
: Particle::Status::Fragmented );
pa.addMother( moth );
}
}
return true;
}
bool
Pythia6Hadroniser::prepareHadronisation( Event& ev )
{
CG_DEBUG( "Pythia6Hadroniser" ) << "Hadronisation preparation called.";
for ( const auto& part : ev.particles() ) {
- //--- loop over all undecayed particles
- if ( part.status() != Particle::Status::Undecayed )
+ if ( part.status() != Particle::Status::Unfragmented )
continue;
+ //--- only loop over all protons to be fragmented
- //--- proton to be fragmented
- const double ranudq = drand();
- short singlet_id = 0, doublet_id = 0;
- if ( ranudq < 1./9. ) {
- singlet_id = PDG::down;
- doublet_id = 2203; // uu1 diquark
- }
- else if ( ranudq < 5./9. ) {
- singlet_id = PDG::up;
- doublet_id = 2101; // ud0 diquark
- }
- else {
- singlet_id = PDG::up;
- doublet_id = 2103; // ud1 diquark
- }
- const double ulmdq = pymass( doublet_id );
- const double ulmq = pymass( singlet_id );
- CG_INFO("") <<ulmdq<<"|"<<ulmq;
+part.dump();
- // Choose random direction in MX frame
- const double ranmxp = 2.*M_PI*drand(); // phi angle
- const double ranmxt = acos( 2.*drand()-1. ); // theta angle
+ const auto partons = pickPartonsContent();
+ const double mx2 = part.mass2();
+ const double mq = pymass( partons.first ), mq2 = mq*mq;
+ const double mdq = pymass( partons.second ), mdq2 = mdq*mdq;
- // Compute momentum of decay particles from MX
- const double pmxp = std::sqrt( std::pow( part.mass2() - ulmdq*ulmdq + ulmq*ulmq, 2 ) / ( 4.*part.mass2() ) - ulmq*ulmq );
+ //--- choose random direction in MX frame
+ const double phi = 2.*M_PI*drand(), theta = acos( 2.*drand()-1. ); // theta angle
- // Build 4-vectors and boost decay particles
+ //--- compute momentum of decay particles from MX
+ const double px = std::sqrt( std::pow( mx2-mdq2+mq2, 2 ) / ( 4.*mx2 ) - mq2 );
- // Start with the singlet
- const auto pmxda_1 = Particle::Momentum::fromPThetaPhi( pmxp, ranmxt, ranmxp, std::hypot( pmxp, ulmq ) ), pmxda_2 = -pmxda_1;
+ //--- build 4-vectors and boost decay particles
+ const auto pq = Particle::Momentum::fromPThetaPhi( px, theta, phi, std::hypot( px, mq ) );
+ //--- singlet
auto singl_mom = part.momentum();
- singl_mom.lorentzBoost( pmxda_1 );
+ singl_mom.lorentzBoost( pq );
- auto& singlet = ev.addParticle( part.role() );
- singlet.setPdgId( singlet_id );
- singlet.setStatus( Particle::Status::DebugResonance );
- singlet.setMomentum( singl_mom );
- //singlet.setMass(); //FIXME
+ auto& quark = ev.addParticle( part.role(), false );
+ quark.setPdgId( partons.first );
+ quark.setStatus( Particle::Status::DebugResonance );
+ quark.setMomentum( singl_mom );
- // Continue with the doublet
+ //--- doublet
auto doubl_mom = part.momentum();
- doubl_mom.lorentzBoost( pmxda_2 );
-
- auto& doublet = ev.addParticle( part.role() );
- doublet.setPdgId( doublet_id );
- doublet.setStatus( Particle::Status::DebugResonance );
- doublet.setMomentum( doubl_mom );
- //std::cout << "doublet, mass = " << doublet.mass() << std::endl;
- //doublet.setMass(); //FIXME
- doublet.dump();
-
- if ( part.numDaughters() == 0 ) {
- singlet.addMother( ev[part.id()] );
- doublet.addMother( ev[part.id()] );
+ doubl_mom.lorentzBoost( -pq );
+
+ auto& diquark = ev.addParticle( part.role(), false );
+ diquark.setPdgId( partons.second );
+ diquark.setStatus( Particle::Status::DebugResonance );
+ diquark.setMomentum( doubl_mom );
+
+ std::cout << (singl_mom+doubl_mom)-part.momentum() << "|" << part.numDaughters() << std::endl;
+ /*if ( part.numDaughters() == 0 ) {
+ quark.addMother( ev[part.id()] );
+ diquark.addMother( ev[part.id()] );
CG_DEBUG( "Pythia6Hadroniser" )
<< "Quark/diquark content succesfully added to the event.";
}
else { // Quark/diquark content already present in the event
CG_WARNING( "Pythia6Hadroniser" )
<< "Quark/diquark content already present in the event!\n\t"
<< "Role of these particles: " << part.role() << ".";
for ( const auto& daug : part.daughters() ) {
if ( ev[daug].pdgId() == PDG::up || ev[daug].pdgId() == PDG::down ) {
//--- quark
- singlet.addMother( ev[part.id()] );
- ev[daug] = singlet;
+ quark.addMother( ev[part.id()] );
+ ev[daug] = quark;
CG_DEBUG( "Pythia6Hadroniser" ) << "Singlet replaced.";
}
else {
//--- diquark
- doublet.addMother( ev[part.id()] );
- ev[daug] = doublet;
+ diquark.addMother( ev[part.id()] );
+ ev[daug] = diquark;
CG_DEBUG( "Pythia6Hadroniser" ) << "Doublet replaced";
}
}
- }
+ }*/
+ ev[part.id()].setStatus( Particle::Status::Fragmented );
}
+ev.dump();
+
return true;
}
Pythia6Hadroniser::EventProperties
Pythia6Hadroniser::fillParticles( const Event& ev ) const
{
pyjets_.n = 0;
//--- initialising the string fragmentation variables
EventProperties out;
out.str_in_evt = 0;
for ( const auto& role : ev.roles() ) { // loop on roles
unsigned int part_in_str = 0;
for ( const auto& part : ev[role] ) {
unsigned int np = part.id();
pyjets_.p[0][np] = part.momentum().px();
pyjets_.p[1][np] = part.momentum().py();
pyjets_.p[2][np] = part.momentum().pz();
pyjets_.p[3][np] = part.energy();
pyjets_.p[4][np] = part.mass();
pyjets_.k[0][np] = part.status() <= Particle::Status::Undefined
? 21 // incoming beam
: (int)part.status();
pyjets_.k[1][np] = part.integerPdgId();
pyjets_.k[2][np] = part.mothers().empty()
? 0 // no mother
: *( part.mothers().begin() )+1; // mother
const auto& daug = part.daughters();
if ( daug.empty() )
pyjets_.k[3][np] = pyjets_.k[4][np] = 0; // no daughters
else {
pyjets_.k[3][np] = *daug.begin()+1; // daughter 1
pyjets_.k[4][np] = *daug.rbegin()+1; // daughter 2
}
for ( int i = 0; i < 5; ++i )
pyjets_.v[i][np] = 0.;
if ( part.status() == Particle::Status::DebugResonance ) {
pyjets_.k[0][np] = 1; //FIXME PYTHIA/JETSET workaround
out.jlpsf[out.str_in_evt][part_in_str] = part.id()+1;
out.num_part_in_str[out.str_in_evt]++;
part_in_str++;
}
pyjets_.n++;
}
if ( out.jlrole[out.str_in_evt] > 0 )
out.str_in_evt++;
}
return out;
}
+
+ std::pair<short,short>
+ Pythia6Hadroniser::pickPartonsContent() const
+ {
+ const double ranudq = drand();
+ if ( ranudq < 1./9. )
+ return { PDG::down, 2203 }; // (d,uu1)
+ if ( ranudq < 5./9. )
+ return { PDG::up, 2101 }; // (u,ud0)
+ return { PDG::up, 2103 }; // (u,ud1)
+ }
+
}
}
// register hadroniser and define alias
REGISTER_HADRONISER( pythia6, Pythia6Hadroniser )
diff --git a/CepGen/Processes/GamGamLL.cpp b/CepGen/Processes/GamGamLL.cpp
index e42fa0c..6061181 100644
--- a/CepGen/Processes/GamGamLL.cpp
+++ b/CepGen/Processes/GamGamLL.cpp
@@ -1,1113 +1,1148 @@
#include "CepGen/Processes/GamGamLL.h"
#include "CepGen/Processes/ProcessesHandler.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"
namespace cepgen
{
namespace proc
{
GamGamLL::GamGamLL( const ParametersList& params ) :
GenericProcess( params, "lpair", "pp → p(*) ( ɣɣ → l⁺l¯ ) p(*)" ),
n_opt_( params.get<int>( "nopt", 0 ) ),
pair_ ( params.get<ParticleProperties>( "pair" ).pdgid ),
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, { pair_, pair_ } }
} );
}
unsigned int
GamGamLL::numDimensions() const
{
switch ( kin_.mode ) {
case KinematicsMode::ElectronProton: default:
throw CG_FATAL( "GamGamLL" )
<< "Process mode " << kin_.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_)[Particle::CentralSystem][0].mass2();
w_limits_ = kin_.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_ = kin_.cuts.initial.q2;
- mx_limits_ = kin_.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;
CG_DEBUG_LOOP( "GamGamLL" )
<< "mc4 = " << mc4_ << "\n\t"
<< "sig1 = " << sig1 << ".";
// 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" );
+ const auto s2 = map( x(2), Limits( sig1, smax ), "s2" );
+ s2_ = s2.first;
+ ds2 = s2.second;
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() );
+ if ( t1_max > -kin_.cuts.initial.q2.min() ) {
+ CG_WARNING( "GamGamLL" ) << "t1max = " << t1_max << " > -q2min = " << -kin_.cuts.initial.q2.min();
return false;
}
- if ( t1_min < -q2_limits_.max() && q2_limits_.hasMax() ) {
- CG_DEBUG( "GamGamLL" ) << "t1min = " << t1_min << " < -q2max = " << -q2_limits_.max();
+ if ( t1_min < -kin_.cuts.initial.q2.max() && kin_.cuts.initial.q2.hasMax() ) {
+ CG_DEBUG( "GamGamLL" ) << "t1min = " << t1_min << " < -q2max = " << -kin_.cuts.initial.q2.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();
+ if ( t1_max < -kin_.cuts.initial.q2.max() && kin_.cuts.initial.q2.hasMax() )
+ t1_max = -kin_.cuts.initial.q2.max();
+ if ( t1_min > -kin_.cuts.initial.q2.min() && kin_.cuts.initial.q2.hasMin() )
+ t1_min = -kin_.cuts.initial.q2.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.;
+ const auto t1 = map( x(0), Limits( t1_min, t1_max ), "t1" );
+ t1_ = t1.first;
+ const double dt1 = -t1.second; // changes wrt mapt1 : dx->-dx
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_ );
Limits s2_lim;
s2_lim.min() = sig*sig;
// one computes splus and (s2x=s2max)
double splus;
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;
s2_lim.max() = se/splus;
}
else {
s2_lim.max() = sb+sd;
splus = se/s2_lim.max();
}
}
else { // 3
s2_lim.max() = ( s_*( t1_*( s_+d8-masses_.MX2 )-w2_*masses_.MX2 )+w2_*masses_.MX2*( w2_+masses_.MX2-t1_ ) )/( ss*t13 );
splus = s2_lim.min();
}
// 4
double s2x = s2_lim.max();
CG_DEBUG_LOOP( "GamGamLL" )
<< "s2x = s2max = " << s2x;
if ( n_opt_ < 0 ) { // 5
if ( splus > s2_lim.min() ) {
s2_lim.min() = splus;
CG_DEBUG_LOOP( "GamGamLL" )
<< "min(sig2) truncated to splus = " << splus;
}
- if ( n_opt_ < -1 )
- map( x(2), s2_lim, s2_, ds2, "s2" );
- else
- mapla( t1_, w2_, x(2), s2_lim, s2_, ds2 ); // n_opt_==-1
+ const auto s2 = n_opt_ < -1
+ ? map( x(2), s2_lim, "s2" )
+ : mapla( t1_, w2_, x(2), s2_lim ); // n_opt_==-1
+ s2_ = s2.first;
+ ds2 = s2.second;
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.;
+ const auto t2 = map( x(1), Limits( t2_min, t2_max ), "t2" );
+ t2_ = t2.first;
+ const double dt2 = -t2.second; // changes wrt mapt2 : dx->-dx
// \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" )
+ << "tau= " << tau << "\n\t"
<< "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";
+ 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;
if ( fabs( ( sl5-b )/sl5 ) >= 1. ) {
s2p = 0.5 * ( sl5-b )/t2_;
s2_lim.min() = c/( t2_*s2p );
}
else { // 8
s2_lim.min() = 0.5 * ( -sl5-b )/t2_;
s2p = c/( t2_*s2_lim.min() );
}
// 9
- if ( n_opt_ > 1 )
- map( x( 2 ), s2_lim, s2_, ds2, "s2" );
- else if ( n_opt_ == 1 )
- mapla( t1_, w2_, x( 2 ), s2_lim, s2_, ds2 );
+ if ( n_opt_ >= 1 ) {
+ const auto s2 = n_opt_ > 1
+ ? map( x( 2 ), s2_lim, "s2" )
+ : mapla( t1_, w2_, x( 2 ), s2_lim );
+ s2_ = s2.first;
+ ds2 = s2.second;
+ }
const double ap = -0.25*pow( s2_+d8, 2 )+s2_*t1_;
dd1_ = 0.25 * ( s2_-s2_lim.max() ) * ( w1_ != 0. ? ( splus-s2_ ) * w1_ : ss * t13 );
dd2_ = 0.25 * ( s2_-s2_lim.min() ) * ( s2p-s2_ ) * t2_;
CG_DEBUG_LOOP( "GamGamLL" )
<< "t2 = " << t2_ << "\n\t"
<< "s2 = " << s2_ << "\n\t"
<< "s2p = " << s2p << "\n\t"
<< "splus = " << splus << "\n\t"
- << "s2 range= " << s2_lim << "\n\t"
- << "dd2 = " << dd2_;;
+ << "s2 range= " << s2_lim;
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 );
+ CG_DEBUG_LOOP( "GamGamLL" )
+ << std::scientific
+ << "dd = " << dd << ", "
+ << "dd1 = " << dd1_ << ", "
+ << "dd2 = " << dd2_ << std::fixed;
+
if ( dd <= 0. ) {
- CG_DEBUG_LOOP( "GamGamLL" )
- << std::scientific
- << "dd = " << dd << " <= 0\n\t"
- << "dd1 = " << dd1_ << "\t"
- << "dd2 = " << dd2_ << std::fixed;
+ CG_WARNING( "GamGamLL:pickin" ) << "dd = " << dd << " <= 0.";
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";
+ CG_WARNING( "GamGamLL:pickin" ) << "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;
+ << "ds2=" << ds2 << ", dt1=" << dt1 << ", dt2=" << dt2 << "\n\t"
+ << "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;
}
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 pp3 = sqrt( ep3*ep3-masses_.MX2 ), pt3 = sqrt( dd1_/s_ )/p_cm_;
const double pp5 = sqrt( ep5*ep5-masses_.MY2 ), pt5 = sqrt( dd3_/s_ )/p_cm_;
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;
}
const double ct3 = ( ep1_*ep3 < p13_ ? -1. : +1. )*sqrt( 1.-sin_theta3*sin_theta3 );
const double ct5 = ( ep2_*ep5 > p25_ ? -1. : +1. )*sqrt( 1.-sin_theta5*sin_theta5 );
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]";
+ << "sin(phi_3) = " << sin_phi3 << " while it must be in (" << Limits( -1., 1. ) << ")";
return false;
}
if ( fabs( sin_phi5 ) > 1. ) {
CG_WARNING( "GamGamLL" )
- << "sin(phi_5) = " << sin_phi5 << " while it must be in [-1 ; 1]";
+ << "sin(phi_5) = " << sin_phi5 << " while it must be in (" << Limits( -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();
+ p5_lab_[0] *= -1.;
else
- p3_lab_[0] = -p3_lab_.px();
+ p3_lab_[0] *= -1.;
return true;
}
//---------------------------------------------------------------------------------------------
double
- GamGamLL::computeOutgoingPrimaryParticlesMasses( double x, double outmass, double lepmass, double& dw )
+ GamGamLL::computeOutgoingPrimaryParticlesMasses( double x, double outmass, double lepmass, double& dmx )
{
- const double mx0 = mp_+PDG::get().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 );
+ const double mx0 = mp_+PDG::get().mass( PDG::piPlus ); // 1.07
+// const double mx0 = 1.07;//FIXME
+ const Limits wx_lim( pow( std::max( mx0, kin_.cuts.remnants.mass_single.min() ), 2 ),
+ pow( std::min( sqs_-outmass-2.*lepmass, kin_.cuts.remnants.mass_single.max() ), 2 ) );
- double mx2 = 0., dmx2 = 0.;
- map( x, Limits( wx2min, wx2max ), mx2, dmx2, "mx2" );
+ const auto w = map( x, wx_lim, "mx2" );
+ const double mx = sqrt( w.first );
+ dmx = sqrt( w.second );
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 );
+ << "mX² in range " << wx_lim << ", x = " << x << "\n\t"
+ << "mX² = " << w.first << ", dmX² = " << w.second << "\n\t"
+ << "mX = " << mx << ", dmX = " << dmx;
- dw = sqrt( dmx2 );
- return sqrt( mx2 );
+ return mx;
}
//---------------------------------------------------------------------------------------------
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();
switch ( kin_.mode ) {
case KinematicsMode::ElectronProton: default:
throw CG_FATAL( "GamGamLL" ) << "Case not yet supported!";
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( PDG::get().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( PDG::get().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" );
+ const auto w4 = map( x( 4 ), w_limits_, "w4" );
+ w4_ = w4.first;
+ const double dw4 = w4.second;
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 = std::hypot( pgx, pgy ), // outgoing proton (3)'s transverse momentum
pgg = std::hypot( pgp, 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;
+ q1dq2_ = 0.5*( w4_-t1_-t2_ );
+
+
+ CG_DEBUG_LOOP( "GamGamLL" )
+ << "ecm6 = " << ecm6 << ", mc4 = " << mc4_ << "\n\t"
+ << "eg = " << eg << ", pg = " << pg << "\n\t"
+ << "q1dq = " << q1dq_ << ", q1dq2 = " << q1dq2_;
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
+ CG_DEBUG_LOOP( "GamGamLL:gmufil" )
+ << "unboosted P(l1)=" << p6_cm_ << "\n\t"
+ << "unboosted P(l2)=" << p7_cm_;
+
p6_cm_.betaGammaBoost( gamma, betgam );
p7_cm_.betaGammaBoost( gamma, betgam );
//--- cut on mass of final hadronic system (MX/Y)
- if ( mx_limits_.valid() ) {
+ if ( kin_.cuts.remnants.mass_single.valid() ) {
if ( ( kin_.mode == KinematicsMode::InelasticElastic
|| kin_.mode == KinematicsMode::InelasticInelastic )
- && !mx_limits_.passes( MX_ ) )
+ && !kin_.cuts.remnants.mass_single.passes( MX_ ) )
return 0.;
if ( ( kin_.mode == KinematicsMode::ElasticInelastic
|| kin_.mode == KinematicsMode::InelasticInelastic )
- && !mx_limits_.passes( MY_ ) )
+ && !kin_.cuts.remnants.mass_single.passes( MY_ ) )
return 0.;
}
//--- cut on the proton's Q2 (first photon propagator T1)
if ( !kin_.cuts.initial.q2.passes( -t1_ ) )
return 0.;
//--- cuts on outgoing leptons' kinematics
if ( !kin_.cuts.central.mass_sum.passes( ( p6_cm_+p7_cm_ ).mass() ) )
return 0.;
//----- cuts on the individual leptons
if ( kin_.cuts.central.pt_single.valid() ) {
const Limits& pt_limits = kin_.cuts.central.pt_single;
if ( !pt_limits.passes( p6_cm_.pt() ) || !pt_limits.passes( p7_cm_.pt() ) )
return 0.;
}
if ( kin_.cuts.central.energy_single.valid() ) {
const Limits& energy_limits = kin_.cuts.central.energy_single;
if ( !energy_limits.passes( p6_cm_.energy() ) || !energy_limits.passes( p7_cm_.energy() ) )
return 0.;
}
if ( kin_.cuts.central.eta_single.valid() ) {
const Limits& eta_limits = kin_.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 ( kin_.mode ) { // inherited from CDF version
case KinematicsMode::ElectronProton: default: jacobian_ *= periPP( 1, 2 ); break;
case KinematicsMode::ElasticElastic: jacobian_ *= periPP( 2, 2 ); break;
case KinematicsMode::InelasticElastic: jacobian_ *= periPP( 3, 2 )*pow( masses_.dw31, 2 ); break;
case KinematicsMode::ElasticInelastic: jacobian_ *= periPP( 3, 2 )*pow( masses_.dw52, 2 ); break;
case KinematicsMode::InelasticInelastic: jacobian_ *= periPP( 3, 3 )*pow( masses_.dw31*masses_.dw52, 2 ); break;
}
+ CG_DEBUG_LOOP( "GamGamLL:f" )
+ << "kinematics mode: " << kin_.mode << "\n\t"
+ << "Jacobian: " << jacobian_;
+
//--- compute the event weight using the Jacobian
return constants::GEVM2_TO_PB*jacobian_;
}
//---------------------------------------------------------------------------------------------
void
GamGamLL::fillKinematics( bool )
{
const Particle::Momentum cm = (*event_)[Particle::IncomingBeam1][0].momentum()
+ (*event_)[Particle::IncomingBeam2][0].momentum();
const double gamma = cm.energy()/sqs_, betgam = cm.pz()/sqs_;
+ CG_DEBUG_LOOP( "GamGamLL:gmufil" )
+ << "sqrt(s)=" << sqs_ << " GeV, initial two-proton system: " << cm << "\n\t"
+ << "gamma=" << gamma << ", betgam=" << betgam;
+
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 );
+
+ CG_DEBUG_LOOP( "GamGamLL:gmufil" )
+ << "unboosted PX=" << p3_lab_ << "\n\t"
+ << "unboosted PY=" << p5_lab_;
+
p3_lab_.betaGammaBoost( gamma, betgam );
p5_lab_.betaGammaBoost( gamma, betgam );
+ CG_DEBUG_LOOP( "GamGamLL:gmufil" )
+ << "boosted PX=" << p3_lab_ << "\n\t"
+ << "boosted PY=" << p5_lab_ << "\n\t"
+ << "boosted P(l1)=" << p6_cm_ << "\n\t"
+ << "boosted P(l2)=" << p7_cm_;
+
//----- 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;
+ const short rany = drand() > 0.5 ? 1 : -1, ransign = drand() > 0.5 ? 1 : -1;
+ const double ranphi = 2*drand()*M_PI;
+
+ Particle::Momentum plab_ph1 = plab_ip1-p3_lab_;
+ plab_ph1.rotatePhi( ranphi, rany );
- Particle::Momentum plab_ph1 = ( plab_ip1-p3_lab_ ).rotatePhi( ranphi, rany );
- Particle::Momentum plab_ph2 = ( plab_ip2-p5_lab_ ).rotatePhi( ranphi, rany );
+ Particle::Momentum plab_ph2 = plab_ip2-p5_lab_;
+ plab_ph2.rotatePhi( ranphi, rany );
p3_lab_.rotatePhi( ranphi, rany );
p5_lab_.rotatePhi( ranphi, rany );
p6_cm_.rotatePhi( ranphi, rany );
p7_cm_.rotatePhi( ranphi, rany );
+ CG_DEBUG_LOOP( "GamGamLL:gmufil" )
+ << "boosted+rotated PX=" << p3_lab_ << "\n\t"
+ << "boosted+rotated PY=" << p5_lab_ << "\n\t"
+ << "boosted+rotated P(l1)=" << p6_cm_ << "\n\t"
+ << "boosted+rotated P(l2)=" << p7_cm_;
+
/*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 );
+ auto& op1 = event_->getOneByRole( Particle::OutgoingBeam1 );
op1.setMomentum( p3_lab_ );
switch ( kin_.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 );
+ auto& op2 = event_->getOneByRole( Particle::OutgoingBeam2 );
op2.setMomentum( p5_lab_ );
switch ( kin_.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 );
+ auto& ph1 = event_->getOneByRole( Particle::Parton1 );
ph1.setMomentum( plab_ph1 );
//----- second incoming photon
- Particle& ph2 = event_->getOneByRole( Particle::Parton2 );
+ auto& ph2 = event_->getOneByRole( Particle::Parton2 );
ph2.setMomentum( plab_ph2 );
- Particles& central_system = (*event_)[Particle::CentralSystem];
+ auto& central_system = (*event_)[Particle::CentralSystem];
//----- first outgoing lepton
- Particle& ol1 = central_system[0];
+ auto& 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];
+ auto& 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" )
+ CG_DEBUG_LOOP( "GamGamLL:peripp" )
<< " Nup = " << nup_ << "\n\t"
<< "Ndown = " << ndown_;
FormFactors fp1, fp2;
formFactors( -t1_, -t2_, fp1, fp2 );
- CG_DEBUG_LOOP( "GamGamLL" )
+ CG_DEBUG_LOOP( "GamGamLL:peripp" )
<< "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" )
+ CG_DEBUG_LOOP( "GamGamLL:peripp" )
+ << "bb = " << bb_ << ", qqq = " << qqq << ", qdq = " << qdq << "\n\t"
<< "t11 = " << t11 << "\t" << "t12 = " << t12 << "\n\t"
<< "t21 = " << t21 << "\t" << "t22 = " << t22 << "\n\t"
- << "⇒ PeriPP = " << peripp;
+ << "=> PeriPP = " << peripp;
return peripp;
}
- void
- GamGamLL::map( double expo, const Limits& lim, double& out, double& dout, const std::string& var_name_ )
+ std::pair<double,double>
+ GamGamLL::map( double expo, const Limits& lim, const std::string& var_name_ )
{
- const double y = lim.max()/lim.min();
- out = lim.min()*pow( y, expo );
- dout = out*log( y );
+ 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;
+ << "Mapping variable \"" << var_name_ << "\" in range (" << lim << ")"
+ << " (max/min = " << y << ")\n\t"
+ << "exponent = " << expo << " => "
+ << "x = " << out << ", dx = " << dout;
+ return { out, dout };
}
- void
- GamGamLL::mapla( double y, double z, int u, const Limits& lim, double& out, double& dout )
+ std::pair<double,double>
+ GamGamLL::mapla( double y, double z, int u, const Limits& lim )
{
const double xmb = lim.min()-y-z, xpb = lim.max()-y-z;
const double c = -4.*y*z;
- const double alp = sqrt( xpb*xpb + c ), alm = sqrt( xmb*xmb + c );
+ 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 );
- out = y+z+( am*zz - c / ( am*zz ) ) / 2.;
+ const double out = y+z+0.5*( am*zz-c / ( am*zz ) );
const double ax = sqrt( pow( out-y-z, 2 )+c );
- dout = ax*log( yy );
+ return { out, ax*log( yy ) };
}
void
GamGamLL::formFactors( double q1, double q2, FormFactors& fp1, FormFactors& fp2 ) const
{
const double mx2 = MX_*MX_, my2 = MY_*MY_;
switch ( kin_.mode ) {
case KinematicsMode::ElasticElastic: default: {
fp1 = FormFactors::protonElastic( -t1_ );
fp2 = FormFactors::protonElastic( -t2_ );
} break;
case KinematicsMode::ElasticInelastic: {
fp1 = FormFactors::protonElastic( -t1_ );
fp2 = FormFactors::protonInelastic( -t2_, w2_, my2, *kin_.structure_functions );
} break;
case KinematicsMode::InelasticElastic: {
fp1 = FormFactors::protonInelastic( -t1_, w1_, mx2, *kin_.structure_functions );
fp2 = FormFactors::protonElastic( -t2_ );
} break;
case KinematicsMode::InelasticInelastic: {
fp1 = FormFactors::protonInelastic( -t1_, w1_, mx2, *kin_.structure_functions );
fp2 = FormFactors::protonInelastic( -t2_, w2_, my2, *kin_.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 d41260a..a0f4555 100644
--- a/CepGen/Processes/GamGamLL.h
+++ b/CepGen/Processes/GamGamLL.h
@@ -1,216 +1,213 @@
#ifndef CepGen_Processes_GamGamLL_h
#define CepGen_Processes_GamGamLL_h
#include "CepGen/Processes/GenericProcess.h"
namespace cepgen
{
namespace proc
{
/**
* 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( *this ) );
}
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
+ /// \param[out] dmx The size of the integration bin
/// \return Mass of the outgoing proton remnant
- double computeOutgoingPrimaryParticlesMasses( double x, double outmass, double lepmass, double& dw );
+ double computeOutgoingPrimaryParticlesMasses( double x, double outmass, double lepmass, double& dmx );
/// 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_;
pdgid_t pair_;
Limits w_limits_;
- Limits q2_limits_;
- Limits mx_limits_;
struct Masses
{
/// squared mass of the first proton-like outgoing particle
double MX2 = 0.;
/// squared mass of the second proton-like outgoing particle
double MY2 = 0.;
/// squared mass of the outgoing leptons
double Ml2 = 0.;
/// \f$\delta_2=m_1^2-m_2^2\f$ as defined in \cite Vermaseren:1982cz
double w12 = 0.;
/// \f$\delta_1=m_3^2-m_1^2\f$ as defined in \cite Vermaseren:1982cz
double w31 = 0.;
double dw31 = 0.;
/// \f$\delta_4=m_5^2-m_2^2\f$ as defined in \cite Vermaseren:1982cz
double w52 = 0.;
double dw52 = 0.;
} 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
+ * \return A pair containing the value and the bin width the new variable definition
* \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, const Limits& lim, double& x, double& d );
+ std::pair<double,double> map( double expo, const Limits& lim, const std::string& var_name = "" );
+ std::pair<double,double> mapla( double y, double z, int u, const Limits& lim );
/// 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/GenericProcess.cpp b/CepGen/Processes/GenericProcess.cpp
index 19e1ae3..155635a 100644
--- a/CepGen/Processes/GenericProcess.cpp
+++ b/CepGen/Processes/GenericProcess.cpp
@@ -1,232 +1,237 @@
#include "CepGen/Processes/GenericProcess.h"
#include "CepGen/Core/Exception.h"
#include "CepGen/Event/Event.h"
#include "CepGen/Physics/Constants.h"
#include "CepGen/Physics/FormFactors.h"
#include "CepGen/Physics/PDG.h"
namespace cepgen
{
namespace proc
{
const double GenericProcess::mp_ = PDG::get().mass( PDG::proton );
const double GenericProcess::mp2_ = GenericProcess::mp_*GenericProcess::mp_;
GenericProcess::GenericProcess( const ParametersList& params, const std::string& name, const std::string& description, bool has_event ) :
params_( params ), 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 ) :
params_( proc.params_ ), 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. ), kin_( proc.kin_ ),
has_event_( proc.has_event_ ), event_( new Event( *proc.event_.get() ) ),
is_point_set_( false )
{}
GenericProcess&
GenericProcess::operator=( const GenericProcess& proc )
{
params_ = proc.params_;
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_;
kin_ = proc.kin_;
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_LOG_MATCH( "Process:dumpPoint", debugInsideLoop ) )
dumpPoint();
clearEvent();
}
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& kin )
{
kin_ = kin;
prepareKinematics();
}
void
GenericProcess::prepareKinematics()
{
if ( !isKinematicsDefined() )
throw CG_FATAL( "GenericProcess" ) << "Kinematics not properly defined for the process.";
const HeavyIon hi1( kin_.incoming_beams.first.pdg ), hi2( kin_.incoming_beams.second.pdg );
const double m1 = hi1 ? HeavyIon::mass( hi1 ) : PDG::get().mass( kin_.incoming_beams.first.pdg );
const double m2 = hi2 ? HeavyIon::mass( hi2 ) : PDG::get().mass( kin_.incoming_beams.second.pdg );
// at some point introduce non head-on colliding beams?
const auto p1 = Particle::Momentum::fromPxPyPzM( 0., 0., +kin_.incoming_beams.first .pz, m1 );
const auto p2 = Particle::Momentum::fromPxPyPzM( 0., 0., -kin_.incoming_beams.second.pz, m2 );
setIncomingKinematics( p1, p2 );
s_ = ( p1+p2 ).mass2();
sqs_ = sqrt( s_ );
w1_ = p1.mass2();
w2_ = p2.mass2();
CG_DEBUG( "GenericProcess" ) << "Kinematics successfully prepared!\n"
<< " √s = " << sqs_*1.e-3 << " TeV,\n"
<< " p₁ = " << p1 << ", mass=" << p1.mass() << " GeV\n"
- << " p₂ = " << p1 << ", mass=" << p2.mass() << " GeV.";
+ << " p₂ = " << p2 << ", mass=" << p2.mass() << " GeV.";
}
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 );
const auto& part_info = PDG::get()( ip.second );
p.setPdgId( ip.second, part_info.charge/3. );
p.setMass( part_info.mass );
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 );
const auto& part_info = PDG::get()( pdg );
p.setPdgId( pdg, part_info.charge/3. );
p.setMass( part_info.mass );
}
}
//----- define the particles parentage
const Particles parts = event_->particles();
for ( const auto& p : parts ) {
Particle& part = (*event_)[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_ )
+ if ( !has_event_ || !event_ )
return;
- event_->getOneByRole( Particle::IncomingBeam1 ).setMomentum( p1 );
- event_->getOneByRole( Particle::IncomingBeam2 ).setMomentum( p2 );
+ CG_DEBUG( "GenericProcess:incomingBeams" )
+ << "Incoming primary particles:\n\t"
+ << p1 << "\n\t"
+ << p2;
+
+ (*event_)[Particle::IncomingBeam1][0].setMomentum( p1 );
+ (*event_)[Particle::IncomingBeam2][0].setMomentum( p2 );
}
bool
GenericProcess::isKinematicsDefined()
{
if ( !has_event_ )
return true;
// check the incoming state
bool is_incoming_state_set =
( !(*event_)[Particle::IncomingBeam1].empty() && !(*event_)[Particle::IncomingBeam2].empty() );
// check the outgoing state
bool is_outgoing_state_set =
( !(*event_)[Particle::OutgoingBeam1].empty() && !(*event_)[Particle::OutgoingBeam2].empty()
&& !(*event_)[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/test/test_singlepoint.cpp b/test/test_singlepoint.cpp
index dbcc63e..ac9f16c 100644
--- a/test/test_singlepoint.cpp
+++ b/test/test_singlepoint.cpp
@@ -1,29 +1,29 @@
#include "CepGen/Cards/Handler.h"
#include "CepGen/Generator.h"
#include "CepGen/Core/Exception.h"
using namespace std;
int main( int argc, char* argv[] )
{
if ( argc < 2 )
throw CG_FATAL( "main" ) << "Usage: " << argv[0] << " input-card";
cepgen::Generator gen;
gen.setParameters( cepgen::card::Handler::parse( argv[1] ) );
CG_INFO( "main" ) << gen.parametersPtr();
cepgen::utils::Logger::get().level = cepgen::utils::Logger::Level::debugInsideLoop;
- vector<double> x( 12, 0.3 );
+ vector<double> x( 12, 0.6 );
cout << "point: ";
string delim;
for ( const auto& v : x )
cout << delim << v, delim = ", ";
cout << endl;
cout << "weight: " << gen.computePoint( &x[0] ) << endl;
return 0;
}