No OneTemporary
Actions

Size

73 KB

Subscribers

None

View Options

	diff --git a/CepGen/Physics/FormFactors.cpp b/CepGen/Physics/FormFactors.cpp
	index 2c7b185..1f0896f 100644
	--- a/CepGen/Physics/FormFactors.cpp
	+++ b/CepGen/Physics/FormFactors.cpp
	@@ -1,109 +1,88 @@
	#include "FormFactors.h"

	namespace CepGen
	{
	FormFactors
	- TrivialFormFactors()
	+ FormFactors::Trivial()
	{
	return FormFactors( 1.0, 1.0 );
	}

	FormFactors
	- ElasticFormFactors( double q2, double mi2 )
	+ FormFactors::ProtonElastic( double q2 )
	{
	- const double GE = pow(1.+q2/0.71, -2.), GM = 2.79*GE;
	- return FormFactors( ( 4.mi2GEGE+q2GMGM ) / ( 4.mi2 + q2 ), GM*GM );
	+ const double mp = Particle::massFromPDGId( Particle::Proton ), mp2 = mp*mp;
	+ const double GE = pow( 1.+q2/0.71, -2. ), GE2 = GE*GE;
	+ const double GM = 2.79GE, GM2 = GMGM;
	+ return FormFactors( ( 4.mp2GE2 + q2GM2 ) / ( 4.mp2 + q2 ), GM2 );
	}

	FormFactors
	- SuriYennieFormFactors( double q2, double mi2, double mf2 )
	+ FormFactors::ProtonInelastic( const StructureFunctionsType& sf, double q2, double mi2, double mf2 )
	{
	+ switch ( sf ) {
	+ case StructureFunctionsType::ElasticProton:
	+ InWarning( "Elastic proton form factors requested! Check your process definition!" );
	+ return FormFactors::ProtonElastic( q2 );
	+ case StructureFunctionsType::SuriYennie:
	+ return FormFactors::SuriYennie( q2, mi2, mf2 );
	+ case StructureFunctionsType::SzczurekUleshchenko:
	+ return FormFactors::SzczurekUleshchenko( q2, mi2, mf2 );
	+ case StructureFunctionsType::Fiore:
	+ case StructureFunctionsType::FioreSea:
	+ case StructureFunctionsType::FioreVal:
	+ return FormFactors::FioreBrasse( q2, mi2, mf2 );
	+ default: throw Exception( __PRETTY_FUNCTION__, "Invalid structure functions required!", FatalError );
	+ }
	+ }
	+
	+ FormFactors
	+ FormFactors::SuriYennie( double q2, double mi2, double mf2 )
	+ {
	+ struct SuriYennieParameters {
	+ SuriYennieParameters( double rho, double bp, double cp, double cc1, double cc2, double dd1 ) :
	+ rho( rho ), Bp( bp ), Cp( cp ), cc1( cc1 ), cc2( cc2 ), dd1( dd1 ) {}
	+ double rho; // in GeV**2
	+ double Bp, Cp;
	+ double cc1, cc2, dd1;
	+ };
	// values extracted from experimental fits
	- const double cc1 = 0.86926, // 0.6303
	- cc2 = 2.23422, // 2.2049
	- dd1 = 0.12549, // 0.0468
	- cp = 0.96, // 1.23
	- bp = 0.63, // 0.61
	- rho = 0.585; // 1.05
	- const double x = q2/(q2+mf2),
	+ SuriYennieParameters sy( 0.585, 0.63, 0.96, 0.86926, 2.23422, 0.12549 );
	+ //SuriYennieParameters sy( 1.05, 0.61, 1.23, 0.6303, 2.2049, 0.0468 );
	+ const double x = q2 / ( q2 + mf2 ),
	dm2 = mf2-mi2,
	- en = dm2+q2,
	- tau = -q2/4./mi2,
	- rhot = rho+q2;
	- const double rho_norm = rho/rhot;
	+ en = q2 + dm2,
	+ tau = q2 / mi2 / 4.,
	+ rhot = sy.rho+q2;
	+ const double rho_norm = sy.rho/rhot;

	FormFactors ff;
	- ff.FM = ( -1./q2 ) * ( -cc1rho_normrho_normdm2 - cc2mi2pow( 1.-x, 4 )/( x( xcp-2bp )+1. ) );
	- ff.FE = ( -tauff.FM + dd1dm2q2rho_normpow( dm2/en, 2 )/( rhotmi2 ) )/( 1. + enen/( 4.mi2*q2 ) );
	+ ff.FM = ( -1./q2 ) * ( -sy.cc1rho_normrho_normdm2 - sy.cc2mi2pow( 1.-x, 4 )/( x( xsy.Cp-2sy.Bp )+1. ) );
	+ ff.FE = ( tauff.FM + sy.dd1dm2q2rho_normpow( dm2/en, 2 )/( rhotmi2 ) )/( 1. + enen/( 4.mi2*q2 ) );
	return ff;
	}

	FormFactors
	- FioreBrasseFormFactors( double q2, double mi2, double mf2 )
	+ FormFactors::FioreBrasse( double q2, double mi2, double mf2 )
	{
	- const double k = 2.*sqrt( mi2 );
	- // start by computing the proton structure function for this Q**2/mX couple
	- double dummy, psfw1, psfw2;
	- if ( !PSF( -q2, mf2, dummy, psfw1, psfw2 ) ) {
	- InWarning( Form( "Fiore-Brasse form factors to be retrieved for an invalid MX value:\n\t"
	- "%.2e GeV, while allowed range is [1.07, 1.99] GeV", sqrt( mf2 ) ) );
	- return FormFactors( 0.0, 0.0 );
	- }
	-
	- return FormFactors( psfw2 / k, -psfw1*k / q2 );
	+ const double x = q2 / ( q2 + mf2 - mi2 ), k = 2.*sqrt( mi2 );
	+ StructureFunctions sf = StructureFunctions::FioreBrasse( q2, x );
	+ return FormFactors( sf.F2 / k, -sf.F1*k / q2 );
	}

	FormFactors
	- SzczurekUleshchenkoFormFactors( double q2, double mi2, double mf2 )
	+ FormFactors::SzczurekUleshchenko( double q2, double mi2, double mf2 )
	{
	- const double k = 2.*sqrt( mi2 ),
	- q2_0 = 0.8;
	-
	- float x = q2 / ( mf2+q2+mi2 ),
	- amu2 = q2+q2_0; // shift the overall scale
	- float xuv, xdv, xus, xds, xss, xg;
	-
	- grv95lo_( x, amu2, xuv, xdv, xus, xds, xss, xg );
	-
	- DebuggingInsideLoop( Form( "Form factor content at xB = %e (scale = %f GeV^2):\n\t"
	- " valence quarks: u / d = %e / %e\n\t"
	- " sea quarks: u / d / s = %e / %e / %e\n\t"
	- " gluons: = %e",
	- x, amu2, xuv, xdv, xus, xds, xss, xg ) );
	-
	- const double F2_aux = 4./9.( xuv+2.xus )
	- + 1./9.( xdv+2.xds )
	- + 1./9.2.xss;
	- /const double F2_aux = 4./9.( xuv + 2.*xus )
	- + 1./9.( 0. + 2.xds )
	- + 1./9.2.xss;*/
	-
	- // F2 corrected for low Q^2 behaviour
	- const double F2_corr = q2/amu2*F2_aux,
	- F1 = F2_corr/( 2.*x ); // Callan-Gross relation
	-
	- /const double F2_corr = Q2 / ( Q2+Q02 ) F2_aux;
	- ///////term1 = pow(1.- x_/2.*(mx2-mp2+Q2)/Q2, 2);
	- //term1 = (1.-x_*(mx2-mp2+Q2)/Q2);
	- term1 = ( 1. - ( Q2-kt2_ ) / Q2 );
	- //term1 = (1.-Q2min/Q2);
	- //term1 = 1.;
	- term2 = pow( kt2_ / ( kt2_+x_(mx2-mp2)+x_x_*mp2 ), 2 );
	- f_aux = F2_corr/( mx2+Q2-mp2 )term1term2;
	- f_ine = Constants::alphaEM/M_PI( 1.-x_ )f_aux/kt2_;
	- return f_ine;*/
	-
	- const double w2 = kx/q2F2_corr,
	- w1 = 2.*F1/k;
	-
	- return FormFactors( w2/k, -w1*k/q2 );
	+ const double x = q2 / ( q2 + mf2 - mi2 );
	+ StructureFunctions sf = StructureFunctions::SzczurekUleshchenko( q2, x );
	+ return FormFactors( sf.F2 * x / q2, -2.*sf.F1 / q2 );
	}

	std::ostream&
	operator<<( std::ostream& os, const FormFactors& ff )
	{
	os << Form( "Form factors: electric: Fe = %.3e ; magnetic: Fm = %.3e", ff.FE, ff.FM ).c_str();
	return os;
	}
	}

	diff --git a/CepGen/Physics/FormFactors.h b/CepGen/Physics/FormFactors.h
	index 5cd6b2c..8b56974 100644
	--- a/CepGen/Physics/FormFactors.h
	+++ b/CepGen/Physics/FormFactors.h
	@@ -1,38 +1,48 @@
	#ifndef CepGen_Physics_FormFactors_h
	#define CepGen_Physics_FormFactors_h

	#include <math.h>
	+#include <array>

	#include "CepGen/Core/utils.h"

	#include "Constants.h"
	#include "StructureFunctions.h"

	namespace CepGen
	{
	/// Form factors collection (electric and magnetic parts)
	- struct FormFactors {
	- /// Initialise a collection of electric/magnetic form factors
	- FormFactors( double fe=0.0, double fm=0.0 ) : FE( fe ), FM( fm ) {}
	- /// Electric form factor
	- double FE;
	- /// Magnetic form factor
	- double FM;
	- /// Dumping operator for standard output streams
	- friend std::ostream& operator<<( std::ostream&, const FormFactors& );
	- };
	+ class FormFactors
	+ {
	+ public:
	+ /// Initialise a collection of electric/magnetic form factors
	+ FormFactors( double fe=0.0, double fm=0.0 ) : FE( fe ), FM( fm ) {}
	+ // compute x from w2/m2
	+ double x( double q2, double w2, double m2=0.0 ) const {
	+ const double mp = Particle::massFromPDGId( Particle::Proton );
	+ return 1./( 1.+( w2-mp*mp ) / q2+m2 );
	+ }
	+ /// Trivial, spin-0 form factors (e.g. pion)
	+ static FormFactors Trivial();
	+ /// Elastic proton form factors
	+ static FormFactors ProtonElastic( double q2 );
	+ /// Suri-Yennie inelastic form factors
	+ static FormFactors SuriYennie( double q2, double mi2, double mf2 );
	+ /// Brasse et al. inelastic form factors
	+ /// \cite Brasse1976413
	+ static FormFactors FioreBrasse( double q2, double mi2, double mf2 );
	+ /// Szczurek-Uleschenko inelastic form factors
	+ static FormFactors SzczurekUleshchenko( double q2, double mi2, double mf2 );
	+ /// Generate the form factors according to the proton structure functions set
	+ static FormFactors ProtonInelastic( const StructureFunctionsType& sf, double q2, double mi2, double mf2 );

	- /// Trivial, spin-0 form factors (e.g. pion)
	- FormFactors TrivialFormFactors();
	- /// Elastic form factors
	- FormFactors ElasticFormFactors( double q2, double mi2 );
	- /// Suri-Yennie inelastic form factors
	- FormFactors SuriYennieFormFactors( double q2, double mi2, double mf2 );
	- /// Brasse et al. inelastic form factors
	- /// \cite Brasse1976413
	- FormFactors FioreBrasseFormFactors( double q2, double mi2, double mf2 );
	- /// Szczurek-Uleschenko inelastic form factors
	- FormFactors SzczurekUleshchenkoFormFactors( double q2, double mi2, double mf2 );
	+ /// Electric form factor
	+ double FE;
	+ /// Magnetic form factor
	+ double FM;
	+ /// Dumping operator for standard output streams
	+ friend std::ostream& operator<<( std::ostream&, const FormFactors& );
	+ };
	}

	#endif
	diff --git a/CepGen/Physics/Kinematics.cpp b/CepGen/Physics/Kinematics.cpp
	index 6fa93db..5c6b375 100644
	--- a/CepGen/Physics/Kinematics.cpp
	+++ b/CepGen/Physics/Kinematics.cpp
	@@ -1,61 +1,59 @@
	#include "Kinematics.h"

	namespace CepGen
	{
	Kinematics::Kinematics() :
	inp( { 6500., 6500. } ), inpdg( { Particle::Proton, Particle::Proton } ),
	central_system( { Particle::Muon, Particle::Muon } ),
	mode( ElasticElastic ), structure_functions( SuriYennie ),
	central_cuts( { { Cuts::pt_single, 3.0 }, { Cuts::pt_diff, { 0., 400.0 } } } ),
	remnant_cuts( { { Cuts::mass, { 1.07, 320.0 } } } ),
	initial_cuts( { { Cuts::q2, { 0.0, 1.0e5 } }, { Cuts::qt, { 0.0, 500.0 } } } )
	{}

	Kinematics::~Kinematics()
	{}

	void
	Kinematics::dump( std::ostream& os ) const
	{
	- os
	- << std::setfill(' ')
	- << __PRETTY_FUNCTION__ << " Dump\n";
	+ os << std::setfill(' ');
	os << "===== Central system\n";
	for ( std::map<Cuts::Central,Limits>::const_iterator lim = central_cuts.begin(); lim != central_cuts.end(); ++lim ) {
	os << std::setw(30) << lim->first << ": " << lim->second;
	}
	os << "===== Initial state\n";
	for ( std::map<Cuts::InitialState,Limits>::const_iterator lim = initial_cuts.begin(); lim != initial_cuts.end(); ++lim ) {
	os << std::setw(30) << lim->first << ": " << lim->second;
	}
	os << "===== Remnants\n";
	for ( std::map<Cuts::Remnants,Limits>::const_iterator lim = remnant_cuts.begin(); lim != remnant_cuts.end(); ++lim ) {
	os << std::setw(30) << lim->first << ": " << lim->second;
	}
	}

	std::ostream&
	operator<<( std::ostream& os, const Kinematics::ProcessMode& pm )
	{
	switch ( pm ) {
	case Kinematics::ElectronElectron: return os << "electron/electron";
	case Kinematics::ElectronProton: return os << "electron/proton";
	case Kinematics::ProtonElectron: return os << "proton/electron";
	case Kinematics::ElasticElastic: return os << "elastic/elastic";
	case Kinematics::InelasticElastic: return os << "inelastic/elastic";
	case Kinematics::ElasticInelastic: return os << "elastic/inelastic";
	case Kinematics::InelasticInelastic: return os << "inelastic/inelastic";
	}
	return os;
	}

	std::ostream&
	operator<<( std::ostream& os, const Kinematics::Limits& lim )
	{
	if ( !lim.hasMin() && !lim.hasMax() ) return os << "no cuts";
	if ( !lim.hasMin() ) return os << Form( "≤ %.3f", lim.max() );
	if ( !lim.hasMax() ) return os << Form( "≥ %.3f", lim.min() );
	return os << Form( "%.3f → %.3f", lim.min(), lim.max() );
	}
	}

	diff --git a/CepGen/Physics/Kinematics.h b/CepGen/Physics/Kinematics.h
	index 107f20b..2ab7ca9 100644
	--- a/CepGen/Physics/Kinematics.h
	+++ b/CepGen/Physics/Kinematics.h
	@@ -1,94 +1,94 @@
	#ifndef CepGen_Physics_Kinematics_h
	#define CepGen_Physics_Kinematics_h

	#include <iomanip>
	#include <algorithm>
	#include <string>

	#include "CepGen/Core/utils.h"
	#include "StructureFunctions.h"
	#include "Particle.h"
	#include "Cuts.h"

	using std::cout;

	namespace CepGen
	{
	/// List of kinematic constraints to apply on the process phase space.
	class Kinematics
	{
	public:
	/// Validity interval for a variable
	class Limits : private std::pair<double,double>
	{
	public:
	/// Define lower and upper limits on a quantity
	Limits( double min=invalid_, double max=invalid_ ) : std::pair<double,double>( min, max ) {}

	/// Lower limit to apply on the variable
	double min() const { return first; }
	/// Lower limit to apply on the variable
	double& min() { return first; }
	/// Upper limit to apply on the variable
	double max() const { return second; }
	/// Upper limit to apply on the variable
	double& max() { return second; }
	/// Specify the lower and upper limits on the variable
	void in( double low, double up ) { first = low; second = up; }
	/// Full variable range allowed
	double range() const { return ( !hasMin() \|\| !hasMax() ) ? 0. : second-first; }
	/// Have a lower limit?
	bool hasMin() const { return first != invalid_; }
	/// Have an upper limit?
	bool hasMax() const { return second != invalid_; }

	/// Human-readable expression of the limits
	friend std::ostream& operator<<( std::ostream&, const Limits& );
	private:
	static constexpr double invalid_ = -999.999;
	};
	public:
	Kinematics();
	~Kinematics();

	/// Type of kinematics to consider for the process
	enum ProcessMode {
	ElectronProton = 0, ///< electron-proton elastic case
	ElasticElastic = 1, ///< proton-proton elastic case
	ElasticInelastic = 2, ///< proton-proton single-dissociative (or inelastic-elastic) case
	InelasticElastic = 3, ///< proton-proton single-dissociative (or elastic-inelastic) case
	InelasticInelastic = 4, ///< proton-proton double-dissociative case
	ProtonElectron,
	ElectronElectron
	};
	/// Human-readable format of a process mode (elastic/dissociative parts)
	friend std::ostream& operator<<( std::ostream&, const ProcessMode& );

	/// Dump all the parameters used in this process cross-section computation
	/// or events generation
	void dump( std::ostream& os=std::cout ) const;

	/// Incoming particles' momentum (in \f$\text{GeV}/c\f$)
	std::pair<double,double> inp;
	/// Set the incoming particles' momenta (if the collision is symmetric)
	inline void setSqrtS( double sqrts ) { inp = { sqrts0.5, sqrts0.5 }; }
	/// Beam/primary particle's PDG identifier
	std::pair<Particle::ParticleCode,Particle::ParticleCode> inpdg;
	/// PDG id of the outgoing central particles
	std::vector<Particle::ParticleCode> central_system;

	/// Type of kinematics to consider for the phase space
	ProcessMode mode;
	/// Type of structure functions to consider
	- StructureFunctions structure_functions;
	+ StructureFunctionsType structure_functions;
	/// Cuts on the central system produced
	std::map<Cuts::Central, Limits> central_cuts;
	/// Cuts on the beam remnants system
	std::map<Cuts::Remnants, Limits> remnant_cuts;
	/// Cuts on the initial particles kinematics
	std::map<Cuts::InitialState, Limits> initial_cuts;
	};
	}

	#endif

	diff --git a/CepGen/Physics/PhotonFluxes.cpp b/CepGen/Physics/PhotonFluxes.cpp
	deleted file mode 100644
	index 92f53e0..0000000
	--- a/CepGen/Physics/PhotonFluxes.cpp
	+++ /dev/null
	@@ -1,73 +0,0 @@
	-#include "PhotonFluxes.h"
	-
	-namespace CepGen
	-{
	- namespace Fluxes
	- {
	- namespace Photon
	- {
	- double
	- ProtonElastic( double x, double kt2 )
	- {
	- const double mp = Particle::massFromPDGId( Particle::Proton ), mp2 = mp*mp;
	-
	- const double Q2_ela = ( kt2+xxmp2 )/( 1.-x );
	- const FormFactors ela = ElasticFormFactors( Q2_ela, mp2 );
	-
	- const double ela1 = pow( kt2/( kt2+xxmp2 ), 2 );
	- const double ela2 = ela.FE;
	- //const double ela3 = 1.-(Q2_ela-kt2)/Q2_ela;
	- //const double ela3 = 1.-xxmp2/Q2_ela/(1.-x);
	- //return Constants::alpha_em/M_PI(1.-x+xx/4.)ela1ela2*ela3/kt2;
	- return Constants::alphaEM/M_PIela1ela2/Q2_ela;
	- //return Constants::alphaEM/M_PI( ( 1.-x )ela1ela2ela3 + xx/2.G_M*G_M )/kt2;
	- }
	-
	-
	- double
	- ProtonInelastic( double x, double kt2, double mx )
	- {
	-#ifndef GRVPDF
	- FatalError( "Inelastic flux cannot be computed as GRV PDF set is not linked to this instance!" );
	-#else
	- const double mx2 = mx*mx,
	- mp = Particle::massFromPDGId( Particle::Proton ),
	- mp2 = mp*mp;
	-
	- const double Q02 = 0.8; // introduced to shift the Q2 scale
	- double term1, term2;
	- double f_aux;
	-
	- // F2 structure function
	- const double Q2min = 1. / ( 1.-x )( x( mx2-mp2 ) + xxmp2 ),
	- Q2 = kt2 / ( 1.-x ) + Q2min;
	- float x_Bjorken = Q2 / ( Q2+mx2-mp2 );
	-
	- float mu2 = Q2+Q02; // scale is shifted
	-
	- float xuv, xdv, xus, xds, xss, xg;
	- grv95lo_( x_Bjorken, mu2, xuv, xdv, xus, xds, xss, xg );
	-
	- const double F2_aux = 4./9.( xuv + 2.xus )
	- + 1./9.( xdv + 2.xds )
	- + 1./9.2.xss;
	-
	- // F2 corrected for low Q^2 behaviour
	- const double F2_corr = Q2 / ( Q2+Q02 ) * F2_aux;
	-
	- ///////term1 = pow(1.- x/2.*(mx2-mp2+Q2)/Q2, 2);
	- //term1 = (1.-x*(mx2-mp2+Q2)/Q2);
	- term1 = ( 1. - ( Q2-kt2 ) / Q2 );
	- //term1 = (1.-Q2min/Q2);
	- //term1 = 1.;
	- term2 = pow( kt2 / ( kt2+x(mx2-mp2)+xx*mp2 ), 2 );
	-
	- f_aux = F2_corr/( mx2+Q2-mp2 )term1term2;
	-
	- return Constants::alphaEM/M_PI( 1.-x )f_aux/kt2;
	-#endif
	- }
	- }
	- }
	-}
	-
	diff --git a/CepGen/Physics/PhotonFluxes.h b/CepGen/Physics/PhotonFluxes.h
	deleted file mode 100644
	index f9cd205..0000000
	--- a/CepGen/Physics/PhotonFluxes.h
	+++ /dev/null
	@@ -1,30 +0,0 @@
	-#ifndef CepGen_Physics_PhotonFluxes_h
	-#define CepGen_Physics_PhotonFluxes_h
	-
	-#include "Particle.h"
	-#include "FormFactors.h"
	-#include "CepGen/Core/Exception.h"
	-
	-namespace CepGen
	-{
	- /// Incoming parton fluxes
	- namespace Fluxes
	- {
	- /// List of fluxes for incoming photons
	- namespace Photon
	- {
	- /// Get the elastic flux to be expected at a given x_bjorken / kT
	- /// \param[in] x Bjorken x
	- /// \param[in] kt2 Transverse 2-momentum \f$\mathbf{q}_{\mathrm{T}}^2\f$ of the incoming photon
	- double ProtonElastic( double x, double kt2 );
	-
	- /// Get the inelastic flux to be expected at a given x_bjorken / kT
	- /// \param[in] x Bjorken x
	- /// \param[in] kt2 Transverse 2-momentum \f$\mathbf{q}_{\mathrm{T}}^2\f$ of the incoming photon
	- /// \param[in] mx Outgoing diffractive proton mass
	- double ProtonInelastic( double x, double kt2, double mx );
	- }
	- }
	-}
	-
	-#endif
	diff --git a/CepGen/Physics/StructureFunctions.cpp b/CepGen/Physics/StructureFunctions.cpp
	index 5298672..aa998d4 100644
	--- a/CepGen/Physics/StructureFunctions.cpp
	+++ b/CepGen/Physics/StructureFunctions.cpp
	@@ -1,92 +1,134 @@
	#include "StructureFunctions.h"

	namespace CepGen
	{
	std::ostream&
	- operator<<( std::ostream& os, const StructureFunctions& sf )
	+ operator<<( std::ostream& os, const StructureFunctionsType& sf )
	{
	switch ( sf ) {
	- case Electron: os << "electron"; break;
	- case ElasticProton: os << "elastic proton"; break;
	- case SuriYennie: os << "Suri-Yennie"; break;
	- case SuriYennieLowQ2: os << "Suri-Yennie;lowQ2"; break;
	- case SzczurekUleshchenko: os << "Szczurek-Uleshchenko"; break;
	- case FioreVal: os << "Fiore;valence"; break;
	- case FioreSea: os << "Fiore;sea"; break;
	- case Fiore: os << "Fiore"; break;
	+ case Electron: return os << "electron";
	+ case ElasticProton: return os << "elastic proton";
	+ case SuriYennie: return os << "Suri-Yennie";
	+ case SuriYennieLowQ2: return os << "Suri-Yennie;lowQ2";
	+ case SzczurekUleshchenko: return os << "Szczurek-Uleshchenko";
	+ case FioreVal: return os << "Fiore;valence";
	+ case FioreSea: return os << "Fiore;sea";
	+ case Fiore: return os << "Fiore";
	+ default: return os;
	}
	- return os;
	}

	/// Fiore-Brasse proton structure function
	- bool
	- PSF( double q2, double mx2, double& sigma_t, double& w1, double& w2 )
	+ StructureFunctions
	+ StructureFunctions::FioreBrasse( double q2, double xbj )
	{
	- sigma_t = w1 = w2 = 0.;
	-
	//const double m_min = Particle::massFromPDGId(Particle::Proton)+0.135;
	- const double m_proton = Particle::massFromPDGId( Particle::Proton ),
	- m2_proton = m_proton*m_proton,
	- m_min = m_proton+Particle::massFromPDGId( Particle::PiZero );
	+ const double mp = Particle::massFromPDGId( Particle::Proton ),mp2 = mp*mp,
	+ m_min = mp+Particle::massFromPDGId( Particle::PiZero );

	- const double mx = sqrt( mx2 );
	+ const double mx2 = mp2 + q2*( 1.-xbj )/xbj, mx = sqrt( mx2 );

	- if ( mx<m_min or mx>1.99 ) {
	- return false;
	+ if ( mx < m_min \|\| mx > 1.99 ) {
	+ InWarning( Form( "Fiore-Brasse form factors to be retrieved for an invalid MX value:\n\t"
	+ "%.2e GeV, while allowed range is [1.07, 1.99] GeV", mx ) );
	+ return StructureFunctions();
	}

	int n_bin;
	double x_bin, dx;
	- if ( mx<1.11 ) {
	+ if ( mx < 1.11 ) {
	n_bin = 0;
	x_bin = mx-m_min;
	dx = 1.11-m_min; // Delta w bin sizes
	}
	- else if ( mx<1.77 ) { // w in [1.11, 1.77[
	+ else if ( mx < 1.77 ) { // w in [1.11, 1.77[
	dx = 0.015; // Delta w bin sizes
	n_bin = ( mx-1.11 )/dx + 1;
	x_bin = fmod( mx-1.11, dx );
	}
	else { // w in [1.77, 1.99[
	dx = 0.02; // Delta w bin sizes
	n_bin = ( mx-1.77 )/dx + 45;
	x_bin = fmod( mx-1.77, dx );
	}

	// values of a, b, c provided from the fits on ep data and retrieved from
	// http://dx.doi.org/10.1016/0550-3213(76)90231-5 with 1.110 <= w2 <=1.990
	- const double abrass[56] = { 5.045, 5.126, 5.390,5.621, 5.913, 5.955,6.139,6.178,6.125, 5.999,
	- 5.769, 5.622, 5.431,5.288, 5.175, 5.131,5.003,5.065,5.045, 5.078,
	- 5.145, 5.156, 5.234,5.298, 5.371, 5.457,5.543,5.519,5.465, 5.384,
	- 5.341, 5.320, 5.275,5.290, 5.330, 5.375,5.428,5.478,5.443, 5.390,
	- 5.333, 5.296, 5.223,5.159, 5.146, 5.143,5.125,5.158,5.159, 5.178,
	- 5.182, 5.195, 5.160,5.195, 5.163, 5.172 },
	- bbrass[56] = { 0.798, 1.052, 1.213,1.334,1.397,1.727,1.750,1.878,1.887,1.927,
	- 2.041, 2.089, 2.148,2.205,2.344,2.324,2.535,2.464,2.564,2.610,
	- 2.609, 2.678, 2.771,2.890,2.982,3.157,3.183,3.315,3.375,3.450,
	- 3.477, 3.471, 3.554,3.633,3.695,3.804,3.900,4.047,4.290,4.519,
	- 4.709, 4.757, 4.840,5.017,5.015,5.129,5.285,5.322,5.545,5.623,
	- 5.775, 5.894, 6.138,6.151,6.301,6.542 },
	- cbrass[56] = { 0.043, 0.024, 0.000,-0.013,-0.023,-0.069,-0.060,-0.080,-0.065,-0.056,
	- -0.065,-0.056,-0.043,-0.034,-0.054,-0.018,-0.046,-0.015,-0.029,-0.048,
	- -0.032,-0.045,-0.084,-0.115,-0.105,-0.159,-0.164,-0.181,-0.203,-0.223,
	- -0.245,-0.254,-0.239,-0.302,-0.299,-0.318,-0.383,-0.393,-0.466,-0.588,
	- -0.622,-0.568,-0.574,-0.727,-0.665,-0.704,-0.856,-0.798,-1.048,-0.980,
	- -1.021,-1.092,-1.313,-1.341,-1.266,-1.473 };
	-
	- const double nu2 = pow( ( mx2-q2-m2_proton ) / ( 2.*m_proton ), 2 ),
	- logqq0 = log( ( nu2-q2 ) / pow( ( mx2-m2_proton ) / ( 2.*m_proton ), 2 ) ) / 2.,
	- gd2 = pow( 1. / ( 1-q2 / .71 ), 4 ); // dipole form factor of the proton
	-
	- const double sigLow = (n_bin==0) ? 0. :
	- exp( abrass[n_bin-1]+bbrass[n_bin-1]logqq0+cbrass[n_bin-1]pow( fabs( logqq0 ), 3 ) )*gd2;
	+ const double a[56] = { 5.045, 5.126, 5.390,5.621, 5.913, 5.955,6.139,6.178,6.125, 5.999,
	+ 5.769, 5.622, 5.431,5.288, 5.175, 5.131,5.003,5.065,5.045, 5.078,
	+ 5.145, 5.156, 5.234,5.298, 5.371, 5.457,5.543,5.519,5.465, 5.384,
	+ 5.341, 5.320, 5.275,5.290, 5.330, 5.375,5.428,5.478,5.443, 5.390,
	+ 5.333, 5.296, 5.223,5.159, 5.146, 5.143,5.125,5.158,5.159, 5.178,
	+ 5.182, 5.195, 5.160,5.195, 5.163, 5.172 },
	+ b[56] = { 0.798, 1.052, 1.213,1.334,1.397,1.727,1.750,1.878,1.887,1.927,
	+ 2.041, 2.089, 2.148,2.205,2.344,2.324,2.535,2.464,2.564,2.610,
	+ 2.609, 2.678, 2.771,2.890,2.982,3.157,3.183,3.315,3.375,3.450,
	+ 3.477, 3.471, 3.554,3.633,3.695,3.804,3.900,4.047,4.290,4.519,
	+ 4.709, 4.757, 4.840,5.017,5.015,5.129,5.285,5.322,5.545,5.623,
	+ 5.775, 5.894, 6.138,6.151,6.301,6.542 },
	+ c[56] = { 0.043, 0.024, 0.000,-0.013,-0.023,-0.069,-0.060,-0.080,-0.065,-0.056,
	+ -0.065,-0.056,-0.043,-0.034,-0.054,-0.018,-0.046,-0.015,-0.029,-0.048,
	+ -0.032,-0.045,-0.084,-0.115,-0.105,-0.159,-0.164,-0.181,-0.203,-0.223,
	+ -0.245,-0.254,-0.239,-0.302,-0.299,-0.318,-0.383,-0.393,-0.466,-0.588,
	+ -0.622,-0.568,-0.574,-0.727,-0.665,-0.704,-0.856,-0.798,-1.048,-0.980,
	+ -1.021,-1.092,-1.313,-1.341,-1.266,-1.473 };
	+
	+ const double d = 3.0;
	+ const double nu = 0.5 * ( q2 + mx2 - mp2 ) / mp, nu2 = nu*nu,
	+ logqq0 = 0.5 * log( ( nu2+q2 ) / pow( ( mx2-mp2 ) / ( 2.*mp ), 2 ) );
	+ const double gd2 = pow( 1. / ( 1+q2 / .71 ), 4 ); // dipole form factor of the proton
	+
	+ const double sigLow = ( n_bin == 0 ) ? 0. :
	+ gd2 * exp( a[n_bin-1] + b[n_bin-1]logqq0 + c[n_bin-1]pow( fabs( logqq0 ), d ) );
	const double sigHigh =
	- exp( abrass[n_bin] +bbrass[n_bin] logqq0+cbrass[n_bin] pow( fabs( logqq0 ), 3 ) )*gd2;
	+ gd2 * exp( a[n_bin] + b[n_bin] logqq0 + c[n_bin] pow( fabs( logqq0 ), d ) );
	+
	+ const double sigma_t = sigLow + x_bin*( sigHigh-sigLow )/dx;
	+ const double w1 = ( mx2-mp2 )/( 8.M_PIM_PImpConstants::alphaEM )/Constants::GeV2toBarn1.e6 sigma_t;
	+ const double w2 = w1 * q2 / ( q2+nu2 );
	+
	+ return StructureFunctions( w1, w2 );
	+ }
	+
	+ StructureFunctions
	+ StructureFunctions::SzczurekUleshchenko( double q2, double xbj )
	+ {
	+#ifndef GRVPDF
	+ FatalError( "Szczurek-Uleshchenko structure functions cannot be computed"
	+ " as GRV PDF set is not linked to this instance!" );
	+#else
	+ const float q02 = 0.8;
	+ float amu2 = q2+q02; // shift the overall scale
	+ float xuv, xdv, xus, xds, xss, xg;
	+ float xbj_arg = xbj;

	- sigma_t = sigLow + x_bin*( sigHigh-sigLow )/dx;
	- w1 = ( mx2-m2_proton )/( 8.M_PIM_PIm_protonConstants::alphaEM )/Constants::GeV2toBarn1.e6 sigma_t;
	- w2 = w1 * q2/( q2-nu2 );
	+ grv95lo_( xbj_arg, amu2, xuv, xdv, xus, xds, xss, xg );

	- return true;
	+ DebuggingInsideLoop( Form( "Form factor content at xB = %e (scale = %f GeV^2):\n\t"
	+ " valence quarks: u / d = %e / %e\n\t"
	+ " sea quarks: u / d / s = %e / %e / %e\n\t"
	+ " gluons: = %e",
	+ xbj, amu2, xuv, xdv, xus, xds, xss, xg ) );
	+
	+ // standard partonic structure function
	+ const double F2_aux = 4./9.( xuv + 2.xus )
	+ + 1./9.( xdv + 2.xds )
	+ + 1./9.( 2.xss );
	+ /const double F2_aux = 4./9.( xuv + 2.*xus )
	+ + 1./9.( 0. + 2.xds )
	+ + 1./9.2.xss;*/
	+
	+ // F2 corrected for low Q^2 behaviour
	+ const double F2_corr = F2_aux * q2 / amu2,
	+ F1 = 0.5*F2_corr/xbj; // Callan-Gross relation
	+
	+ return StructureFunctions( F1, F2_corr );
	+#endif
	+ }
	+
	+ std::ostream&
	+ operator<<( std::ostream& os, const StructureFunctions& sf )
	+ {
	+ return os << "F1 = " << sf.F1 << ", F2 = " << sf.F2;
	}
	}
	diff --git a/CepGen/Physics/StructureFunctions.h b/CepGen/Physics/StructureFunctions.h
	index ec9ce67..582341b 100644
	--- a/CepGen/Physics/StructureFunctions.h
	+++ b/CepGen/Physics/StructureFunctions.h
	@@ -1,41 +1,45 @@
	#ifndef CepGen_Physics_StructureFunctions_h
	#define CepGen_Physics_StructureFunctions_h

	#include "Particle.h"

	extern "C"
	{
	extern void grv95lo_( float&, float&, float&, float&, float&, float&, float&, float& );
	}

	namespace CepGen
	{
	/// Proton structure function to be used in the outgoing state description
	- enum StructureFunctions {
	+ /// \note Values correspond to the LPAIR legacy steering card values
	+ enum StructureFunctionsType {
	Electron = 1,
	ElasticProton = 2,
	SuriYennie = 11,
	SuriYennieLowQ2 = 12,
	SzczurekUleshchenko = 15,
	FioreVal = 101,
	FioreSea = 102,
	Fiore = 103
	};
	+ /// Human-readable format of a structure function type
	+ std::ostream& operator<<( std::ostream& os, const StructureFunctionsType& sf );
	+
	+ class StructureFunctions
	+ {
	+ public:
	+ StructureFunctions( double f1=0.0, double f2=0.0 ) : F1( f1 ), F2( f2 ) {}
	+ /// Fiore-Brasse proton structure functions (F.W Brasse et al., DESY 76/11 (1976),
	+ /// http://dx.doi.org/10.1016/0550-3213(76)90231-5)
	+ /// \param[in] q2 Squared 4-momentum transfer
	+ /// \param[in] xbj Bjorken's x
	+ /// \cite Brasse1976413
	+ static StructureFunctions FioreBrasse( double q2, double xbj );
	+ static StructureFunctions SzczurekUleshchenko( double q2, double xbj );
	+ double F1, F2;
	+ };
	/// Human-readable format of a structure function object
	std::ostream& operator<<( std::ostream& os, const StructureFunctions& sf );
	-
	- /**
	- * Compute the proton structure function (F.W Brasse et al., DESY 76/11 (1976),
	- * http://dx.doi.org/10.1016/0550-3213(76)90231-5)
	- * \param[in] q2 Squared 4-momentum transfer
	- * \param[in] mx2 Squared mass of the proton remnant
	- * \param[out] sigma_t ...
	- * \param[out] w1 First proton structure function: \f$\mathcal W_1\f$
	- * \param[out] w2 Second proton structure function: \f$\mathcal W_2\f$
	- * \cite Brasse1976413
	- */
	- bool PSF( double q2, double mx2, double& sigma_t, double& w1, double& w2 );
	-
	}

	#endif
	diff --git a/CepGen/Processes/GenericKTProcess.cpp b/CepGen/Processes/GenericKTProcess.cpp
	index 4e5cdd6..a43652e 100644
	--- a/CepGen/Processes/GenericKTProcess.cpp
	+++ b/CepGen/Processes/GenericKTProcess.cpp
	@@ -1,227 +1,266 @@
	#include "GenericKTProcess.h"

	namespace CepGen
	{
	namespace Process
	{
	GenericKTProcess::GenericKTProcess( const std::string& name,
	const std::string& description,
	const unsigned int& num_user_dimensions,
	const Particle::ParticleCode& parton1,
	const Particle::ParticleCode& central1,
	const Particle::ParticleCode& parton2,
	const Particle::ParticleCode& central2 ) :
	GenericProcess( name, description+" (kT-factorisation approach)" ),
	kNumUserDimensions( num_user_dimensions ),
	kIntermediatePart1( parton1 ), kIntermediatePart2( parton2 ),
	kProducedPart1( central1 ), kProducedPart2( central2 )
	{
	if ( parton2 == Particle::invalidParticle ) kIntermediatePart2 = kIntermediatePart1;
	if ( central2 == Particle::invalidParticle ) kProducedPart2 = kProducedPart1;
	}

	GenericKTProcess::~GenericKTProcess()
	{}

	void
	GenericKTProcess::addEventContent()
	{
	GenericProcess::setEventContent(
	{ // incoming state
	{ Particle::IncomingBeam1, Particle::Proton },
	{ Particle::IncomingBeam2, Particle::Proton },
	{ Particle::Parton1, kIntermediatePart1 },
	{ Particle::Parton2, kIntermediatePart2 }
	},
	{ // outgoing state
	{ Particle::OutgoingBeam1, { Particle::Proton } },
	{ Particle::OutgoingBeam2, { Particle::Proton } },
	{ Particle::CentralSystem, { kProducedPart1, kProducedPart2 } }
	}
	);
	}

	unsigned int
	GenericKTProcess::numDimensions( const Kinematics::ProcessMode& process_mode ) const
	{
	switch ( process_mode ) {
	default:
	case Kinematics::ElasticElastic: return kNumRequiredDimensions+kNumUserDimensions;
	case Kinematics::ElasticInelastic:
	case Kinematics::InelasticElastic: return kNumRequiredDimensions+kNumUserDimensions+1;
	case Kinematics::InelasticInelastic: return kNumRequiredDimensions+kNumUserDimensions+2;
	}
	}

	void
	GenericKTProcess::addPartonContent()
	{
	// Incoming partons
	qt1_ = exp( log_qmin_+( log_qmax_-log_qmin_ )*x( 0 ) );
	qt2_ = exp( log_qmin_+( log_qmax_-log_qmin_ )*x( 1 ) );
	phi_qt1_ = 2.M_PIx( 2 );
	phi_qt2_ = 2.M_PIx( 3 );
	DebuggingInsideLoop( Form( "photons transverse virtualities (qt):\n\t"
	" mag = %f / %f (%.2f < log(qt) < %.2f)\n\t"
	" phi = %f / %f",
	qt1_, qt2_, log_qmin_, log_qmax_, phi_qt1_, phi_qt2_ ) );
	}

	void
	GenericKTProcess::setKinematics( const Kinematics& kin )
	{
	cuts_ = kin;
	log_qmin_ = -10.; // FIXME //log_qmin_ = std::log( std::sqrt( cuts_.q2min ) );
	log_qmax_ = log( cuts_.initial_cuts[Cuts::qt].max() );
	}

	double
	GenericKTProcess::computeWeight()
	{
	addPartonContent();
	prepareKTKinematics();
	computeOutgoingPrimaryParticlesMasses();

	const double jac = computeJacobian(),
	integrand = computeKTFactorisedMatrixElement(),
	weight = jac*integrand;
	DebuggingInsideLoop( Form( "Jacobian = %f\n\tIntegrand = %f\n\tdW = %f", jac, integrand, weight ) );

	return weight;
	}

	void
	GenericKTProcess::computeOutgoingPrimaryParticlesMasses()
	{
	const unsigned int op_index = kNumRequiredDimensions+kNumUserDimensions;
	const Kinematics::Limits remn_mx_cuts = cuts_.remnant_cuts[Cuts::mass];
	switch ( cuts_.mode ) {
	case Kinematics::ElectronProton: default: {
	InError( "This kT factorisation process is intended for p-on-p collisions! Aborting!" );
	exit( 0 ); } break;
	case Kinematics::ElasticElastic: {
	MX_ = event_->getOneByRole( Particle::IncomingBeam1 ).mass();
	MY_ = event_->getOneByRole( Particle::IncomingBeam2 ).mass();
	} break;
	case Kinematics::ElasticInelastic: {
	const double mx_min = remn_mx_cuts.min(), mx_range = remn_mx_cuts.range();
	MX_ = event_->getOneByRole( Particle::IncomingBeam1 ).mass();
	MY_ = mx_min + mx_range*x( op_index );
	} break;
	case Kinematics::InelasticElastic: {
	const double mx_min = remn_mx_cuts.min(), mx_range = remn_mx_cuts.range();
	MX_ = mx_min + mx_range*x( op_index );
	MY_ = event_->getOneByRole( Particle::IncomingBeam2 ).mass();
	} break;
	case Kinematics::InelasticInelastic: {
	const double mx_min = remn_mx_cuts.min(), mx_range = remn_mx_cuts.range();
	MX_ = mx_min + mx_range*x( op_index );
	MY_ = mx_min + mx_range*x( op_index+1 );
	} break;
	}
	DebuggingInsideLoop( Form( "outgoing remnants invariant mass: %f / %f (%.2f < M(X/Y) < %.2f)", MX_, MY_, remn_mx_cuts.min(), remn_mx_cuts.max() ) );
	}

	void
	GenericKTProcess::computeIncomingFluxes( double x1, double q1t2, double x2, double q2t2 )
	{
	flux1_ = flux2_ = 0.;
	switch ( cuts_.mode ) {
	case Kinematics::ElasticElastic:
	- flux1_ = Fluxes::Photon::ProtonElastic( x1, q1t2 );
	- flux2_ = Fluxes::Photon::ProtonElastic( x2, q2t2 );
	+ flux1_ = elasticFlux( x1, q1t2 );
	+ flux2_ = elasticFlux( x2, q2t2 );
	break;
	case Kinematics::ElasticInelastic:
	- flux1_ = Fluxes::Photon::ProtonElastic( x1, q1t2 );
	- flux2_ = Fluxes::Photon::ProtonInelastic( x2, q2t2, MY_ );
	+ flux1_ = elasticFlux( x1, q1t2 );
	+ flux2_ = inelasticFlux( x2, q2t2, MY_ );
	break;
	case Kinematics::InelasticElastic:
	- flux1_ = Fluxes::Photon::ProtonInelastic( x1, q1t2, MX_ );
	- flux2_ = Fluxes::Photon::ProtonElastic( x2, q2t2 );
	+ flux1_ = inelasticFlux( x1, q1t2, MX_ );
	+ flux2_ = elasticFlux( x2, q2t2 );
	break;
	case Kinematics::InelasticInelastic:
	- flux1_ = Fluxes::Photon::ProtonInelastic( x1, q1t2, MX_ );
	- flux2_ = Fluxes::Photon::ProtonInelastic( x2, q2t2, MY_ );
	+ flux1_ = inelasticFlux( x1, q1t2, MX_ );
	+ flux2_ = inelasticFlux( x2, q2t2, MY_ );
	break;
	default: return;
	}
	if ( flux1_<1.e-20 ) flux1_ = 0.;
	if ( flux2_<1.e-20 ) flux2_ = 0.;
	DebuggingInsideLoop( Form( "Form factors: %e / %e", flux1_, flux2_ ) );
	}

	void
	GenericKTProcess::fillKinematics( bool )
	{
	fillPrimaryParticlesKinematics();
	fillCentralParticlesKinematics();
	}

	void
	GenericKTProcess::fillPrimaryParticlesKinematics()
	{
	//=================================================================
	// outgoing protons
	//=================================================================
	Particle& op1 = event_->getOneByRole( Particle::OutgoingBeam1 ),
	&op2 = event_->getOneByRole( Particle::OutgoingBeam2 );
	// off-shell particles (remnants?)
	bool os1 = false, os2 = false;
	switch ( cuts_.mode ) {
	case Kinematics::ElectronProton: default: {
	InError( "This kT factorisation process is intended for p-on-p collisions! Aborting!" );
	exit(0); } break;
	case Kinematics::ElasticElastic:
	op1.setStatus( Particle::FinalState );
	op2.setStatus( Particle::FinalState );
	break;
	case Kinematics::ElasticInelastic:
	op1.setStatus( Particle::FinalState );
	op2.setStatus( Particle::Undecayed ); op2.setMass(); os2 = true;
	break;
	case Kinematics::InelasticElastic:
	op1.setStatus( Particle::Undecayed ); op1.setMass(); os1 = true;
	op2.setStatus( Particle::FinalState );
	break;
	case Kinematics::InelasticInelastic:
	op1.setStatus( Particle::Undecayed ); op1.setMass(); os1 = true;
	op2.setStatus( Particle::Undecayed ); op2.setMass(); os2 = true;
	break;
	}

	op1.setMomentum( PX_, os1 );
	op2.setMomentum( PY_, os2 );

	//=================================================================
	// incoming partons (photons, pomerons, ...)
	//=================================================================
	//FIXME ensure the validity of this approach
	Particle& g1 = event_->getOneByRole( Particle::Parton1 ),
	&g2 = event_->getOneByRole( Particle::Parton2 );
	g1.setMomentum( event_->getOneByRole( Particle::IncomingBeam1 ).momentum()-PX_, true);
	g1.setStatus( Particle::Incoming );
	g2.setMomentum( event_->getOneByRole( Particle::IncomingBeam2 ).momentum()-PY_, true);
	g2.setStatus( Particle::Incoming );
	}

	double
	GenericKTProcess::minimalJacobian() const
	{
	double jac = 1.;
	jac = ( log_qmax_-log_qmin_ )qt1_; // d(q1t) . q1t
	jac = ( log_qmax_-log_qmin_ )qt2_; // d(q2t) . q2t
	jac = 2.M_PI; // d(phi1)
	jac = 2.M_PI; // d(phi2)

	const double mx_range = cuts_.remnant_cuts.at( Cuts::mass ).range();
	switch ( cuts_.mode ) {
	case Kinematics::ElasticElastic: default: break;
	case Kinematics::ElasticInelastic: jac = 2. mx_range * MY_; break;
	case Kinematics::InelasticElastic: jac = 2. mx_range * MX_; break;
	case Kinematics::InelasticInelastic: jac = 2. mx_range * MX_;
	jac = 2. mx_range * MY_; break;
	} // d(mx/y**2)
	return jac;
	}
	+
	+ double
	+ GenericKTProcess::elasticFlux( double x, double kt2 )
	+ {
	+ const double mp = Particle::massFromPDGId( Particle::Proton ), mp2 = mp*mp;
	+
	+ const double Q2_ela = ( kt2+xxmp2 )/( 1.-x );
	+ const FormFactors ela = FormFactors::ProtonElastic( Q2_ela );
	+
	+ const double ela1 = pow( kt2/( kt2+xxmp2 ), 2 );
	+ const double ela2 = ela.FE;
	+ //const double ela3 = 1.-(Q2_ela-kt2)/Q2_ela;
	+ //const double ela3 = 1.-xxmp2/Q2_ela/(1.-x);
	+ //return Constants::alpha_em/M_PI(1.-x+xx/4.)ela1ela2*ela3/kt2;
	+ return Constants::alphaEM/M_PIela1ela2/Q2_ela;
	+ //return Constants::alphaEM/M_PI( ( 1.-x )ela1ela2ela3 + xx/2.G_M*G_M )/kt2;
	+ }
	+
	+
	+ double
	+ GenericKTProcess::inelasticFlux( double x, double kt2, double mx )
	+ {
	+ const double mx2 = mx*mx,
	+ mp = Particle::massFromPDGId( Particle::Proton ), mp2 = mp*mp;
	+
	+ // F2 structure function
	+ const double Q2min = 1. / ( 1.-x )( x( mx2-mp2 ) + xxmp2 ),
	+ Q2 = kt2 / ( 1.-x ) + Q2min;
	+ float x_Bjorken = Q2 / ( Q2+mx2-mp2 );
	+
	+ StructureFunctions sf = StructureFunctions::SzczurekUleshchenko( Q2, x_Bjorken );
	+
	+ const double term1 = ( 1. - ( Q2-kt2 ) / Q2 ),
	+ term2 = pow( kt2 / ( kt2+x(mx2-mp2)+xx*mp2 ), 2 );
	+
	+ const double f_aux = sf.F2/( mx2+Q2-mp2 )term1term2;
	+
	+ return Constants::alphaEM/M_PI( 1.-x )f_aux/kt2;
	+ }
	}
	}
	diff --git a/CepGen/Processes/GenericKTProcess.h b/CepGen/Processes/GenericKTProcess.h
	index 23ea8da..1406645 100644
	--- a/CepGen/Processes/GenericKTProcess.h
	+++ b/CepGen/Processes/GenericKTProcess.h
	@@ -1,115 +1,126 @@
	#ifndef CepGen_Processes_GenericKTProcess_h
	#define CepGen_Processes_GenericKTProcess_h

	#include "GenericProcess.h"
	#include "CepGen/Physics/FormFactors.h"
	-#include "CepGen/Physics/PhotonFluxes.h"

	namespace CepGen
	{
	namespace Process
	{
	/**
	* A generic kT-factorisation process.
	* First 4 dimensions of the phase space are required for the incoming partons' virtualities (radial and azimuthal coordinates)
	* \brief Class template to define any kT-factorisation process
	* \author Laurent Forthomme <laurent.forthomme@cern.ch>
	* \date Apr 2016
	*/
	class GenericKTProcess : public GenericProcess
	{
	public:
	/**
	* \brief Class constructor
	* \param[in] name Generic process name
	* \param[in] description Human-readable kT-factorised process name
	* \param[in] num_user_dimensions Number of additional dimensions required for the user process
	* \param[in] ip1 First incoming parton
	* \param[in] ip2 Second incoming parton (if undefined, same as the first)
	* \param[in] op1 First produced final state particle
	* \param[in] op2 Second produced final state particle (if undefined, same as the first)
	*/
	GenericKTProcess( const std::string& name,
	const std::string& description="<generic process>",
	const unsigned int& num_user_dimensions=0,
	const Particle::ParticleCode& ip1=Particle::Photon,
	const Particle::ParticleCode& op1=Particle::Muon,
	const Particle::ParticleCode& ip2=Particle::invalidParticle,
	const Particle::ParticleCode& op2=Particle::invalidParticle);
	~GenericKTProcess();

	/// Populate the event content with the generated process' topology
	void addEventContent();
	/// Retrieve the total number of dimensions on which the integration is being performet
	/// \param[in] proc_mode_ Kinematics case considered
	unsigned int numDimensions( const Kinematics::ProcessMode& proc_mode_ ) const;
	/// Retrieve the event weight in the phase space
	double computeWeight();
	/// Populate the event content with the generated process' kinematics
	void fillKinematics( bool );

	protected:
	/// Set the kinematics associated to the phase space definition
	void setKinematics( const Kinematics& kin );
	/// Set the kinematics of the central system before any point computation
	inline virtual void prepareKTKinematics() { DebuggingInsideLoop("Dummy kinematics prepared!"); }
	/// Minimal Jacobian weight of the point considering a kT factorisation
	double minimalJacobian() const;
	/// Jacobian weight of the point in the phase space for integration
	virtual double computeJacobian() = 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 invariant masses of the outgoing protons (or remnants)
	void computeOutgoingPrimaryParticlesMasses();
	/// Compute the unintegrated photon fluxes (for inelastic distributions, interpolation on double logarithmic grid)
	void computeIncomingFluxes( double, double, double, double );
	/// 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;

	+
	+ /// Get the elastic flux to be expected at a given x_bjorken / kT
	+ /// \param[in] x Bjorken x
	+ /// \param[in] kt2 Transverse 2-momentum \f$\mathbf{q}_{\mathrm{T}}^2\f$ of the incoming photon
	+ double elasticFlux( double x, double kt2 );
	+
	+ /// Get the inelastic flux to be expected at a given x_bjorken / kT
	+ /// \param[in] x Bjorken x
	+ /// \param[in] kt2 Transverse 2-momentum \f$\mathbf{q}_{\mathrm{T}}^2\f$ of the incoming photon
	+ /// \param[in] mx Outgoing diffractive proton mass
	+ double inelasticFlux( double x, double kt2, double mx );
	+
	/// Retrieve a component of the phase space point for the kT-factorised process
	inline double xkt( const unsigned int i ) const { return x( kNumRequiredDimensions + i ); }

	/// Minimal log-virtuality of the intermediate parton
	double log_qmin_;
	/// Maximal log-virtuality of the intermediate parton
	double log_qmax_;
	/// 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_;
	/// First incoming parton's flux
	double flux1_;
	/// Second incoming parton's flux
	double flux2_;

	private:
	void addPartonContent();
	static const unsigned int kNumRequiredDimensions = 4;
	/// Number of additional dimensions required for the user process
	/// (in addition to the 4 required for the two partons' transverse momenta)
	unsigned int kNumUserDimensions;
	/// First intermediate parton (photon, pomeron, ...)
	Particle::ParticleCode kIntermediatePart1;
	/// Second intermediate parton (photon, pomeron, ...)
	Particle::ParticleCode kIntermediatePart2;
	/// Type of particle produced in the final state
	Particle::ParticleCode kProducedPart1;
	/// Type of particle produced in the final state
	Particle::ParticleCode kProducedPart2;
	};
	}
	}

	#endif
	diff --git a/CepGen/Processes/GenericProcess.cpp b/CepGen/Processes/GenericProcess.cpp
	index dc82468..f0a1c3a 100644
	--- a/CepGen/Processes/GenericProcess.cpp
	+++ b/CepGen/Processes/GenericProcess.cpp
	@@ -1,202 +1,184 @@
	#include "GenericProcess.h"

	namespace CepGen
	{
	namespace Process
	{
	GenericProcess::GenericProcess( const std::string& name, const std::string& description, bool has_event ) :
	event_( std::shared_ptr<Event>( new Event ) ),
	is_point_set_( false ), is_incoming_state_set_( false ), is_outgoing_state_set_( false ), is_kinematics_set_( false ),
	name_( name ), description_( description ),
	total_gen_time_( 0. ), num_gen_events_( 0 ), has_event_( has_event )
	{}

	GenericProcess::~GenericProcess()
	{}

	void
	GenericProcess::setPoint( const unsigned int ndim, double* x )
	{
	if ( ndim != x_.size() ) x_.resize( ndim );

	x_ = std::vector<double>( x, x+ndim );
	is_point_set_ = true;
	if ( Logger::get().level>=Logger::DebugInsideLoop ) { dumpPoint( DebugMessage ); }
	}

	void
	GenericProcess::prepareKinematics()
	{
	if ( !isKinematicsDefined() ) return; // FIXME dump some information...
	const Particle ib1 = event_->getOneByRole( Particle::IncomingBeam1 ),
	ib2 = event_->getOneByRole( Particle::IncomingBeam2 );

	sqs_ = CMEnergy( ib1, ib2 );
	s_ = sqs_*sqs_;

	w1_ = ib1.mass2();
	w2_ = ib2.mass2();

	Debugging( Form( "Kinematics successfully prepared! sqrt(s) = %.2f", sqs_ ) );
	}

	void
	GenericProcess::dumpPoint( const ExceptionType& et )
	{
	std::ostringstream os;
	for ( unsigned int i=0; i<x_.size(); i++ ) {
	os << Form( " x(%2d) = %8.6f\n\t", i, x_[i] );
	}
	if ( et < DebugMessage ) {
	Information( Form( "Number of integration parameters: %d\n\t"
	"%s", x_.size(), os.str().c_str() ) ); }
	else {
	Debugging( Form( "Number of integration parameters: %d\n\t"
	"%s", x_.size(), os.str().c_str() ) ); }
	}

	void
	GenericProcess::setEventContent( const IncomingState& is, const OutgoingState& os )
	{
	event_->clear();
	//----- add the particles in the event

	//--- incoming state
	for ( IncomingState::const_iterator ip=is.begin(); ip!=is.end(); ip++ ) {
	event_->addParticle( Particle( ip->first, ip->second ) );
	}
	//--- central system (if not already there)
	IncomingState::const_iterator central_system = is.find( Particle::CentralSystem );
	if ( central_system == is.end() ) {
	event_->addParticle( Particle( Particle::Intermediate, Particle::invalidParticle, Particle::Propagator ) );
	}
	//--- outgoing state
	for ( OutgoingState::const_iterator op = os.begin(); op != os.end(); ++op ) {
	for ( std::vector<Particle::ParticleCode>::const_iterator it = op->second.begin(); it != op->second.end(); ++it ) {
	event_->addParticle( Particle( op->first, *it ) );
	}
	}

	//----- define the particles parentage

	const Particles parts = event_->particles();
	for ( Particles::const_iterator p = parts.begin(); p != parts.end(); ++p ) {
	Particle& part = event_->getById( 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_->init();
	}

	void
	GenericProcess::setIncomingKinematics( const Particle::Momentum& p1, const Particle::Momentum& p2 )
	{
	event_->getOneByRole( Particle::IncomingBeam1 ).setMomentum( p1 );
	event_->getOneByRole( Particle::IncomingBeam2 ).setMomentum( p2 );
	}

	void
	GenericProcess::formFactors( double q1, double q2, FormFactors& fp1, FormFactors& fp2 ) const
	{
	const double mx2 = MX_MX_, my2 = MY_MY_;

	- bool inel_p1 = false, inel_p2 = false;
	-
	switch ( cuts_.mode ) {
	case Kinematics::ElectronElectron: {
	- fp1 = TrivialFormFactors(); // electron (trivial) form factor
	- fp2 = TrivialFormFactors(); // electron (trivial) form factor
	+ fp1 = FormFactors::Trivial(); // electron (trivial) form factor
	+ fp2 = FormFactors::Trivial(); // electron (trivial) form factor
	} break;
	case Kinematics::ProtonElectron: {
	- fp1 = ElasticFormFactors( -t1_, w1_ ); // proton elastic form factor
	- fp2 = TrivialFormFactors(); // electron (trivial) form factor
	+ fp1 = FormFactors::ProtonElastic( -t1_ ); // proton elastic form factor
	+ fp2 = FormFactors::Trivial(); // electron (trivial) form factor
	} break;
	case Kinematics::ElectronProton: {
	- fp1 = TrivialFormFactors(); // electron (trivial) form factor
	- fp2 = ElasticFormFactors( -t2_, w2_ ); // proton elastic form factor
	+ fp1 = FormFactors::Trivial(); // electron (trivial) form factor
	+ fp2 = FormFactors::ProtonElastic( -t2_ ); // proton elastic form factor
	} break;
	case Kinematics::ElasticElastic: {
	- fp1 = ElasticFormFactors( -t1_, w1_ ); // proton elastic form factor
	- fp2 = ElasticFormFactors( -t2_, w2_ ); // proton elastic form factor
	+ fp1 = FormFactors::ProtonElastic( -t1_ ); // proton elastic form factor
	+ fp2 = FormFactors::ProtonElastic( -t2_ ); // proton elastic form factor
	} break;
	case Kinematics::ElasticInelastic: {
	- fp1 = ElasticFormFactors( -t1_, w1_ );
	- inel_p2 = true;
	+ fp1 = FormFactors::ProtonElastic( -t1_ );
	+ fp2 = FormFactors::ProtonInelastic( cuts_.structure_functions, -t2_, w2_, my2 );
	} break;
	case Kinematics::InelasticElastic: {
	- inel_p1 = true;
	- fp2 = ElasticFormFactors( -t2_, w2_ );
	+ fp1 = FormFactors::ProtonInelastic( cuts_.structure_functions, -t1_, w1_, mx2 );
	+ fp2 = FormFactors::ProtonElastic( -t2_ );
	} break;
	case Kinematics::InelasticInelastic: {
	- inel_p1 = inel_p2 = true;
	- } break;
	- }
	- switch ( cuts_.structure_functions ) {
	- case SuriYennie:
	- default: {
	- if ( inel_p1 ) fp1 = SuriYennieFormFactors( -t1_, w1_, mx2 );
	- if ( inel_p2 ) fp2 = SuriYennieFormFactors( -t2_, w2_, my2 );
	- } break;
	- case Fiore:
	- case FioreSea:
	- case FioreVal: { // low-Q2 inelastic form factor
	- if ( inel_p1 ) fp1 = FioreBrasseFormFactors( -t1_, w1_, mx2 );
	- if ( inel_p2 ) fp2 = FioreBrasseFormFactors( -t2_, w2_, my2 );
	- } break;
	- case SzczurekUleshchenko: {
	- if ( inel_p1 ) fp1 = SzczurekUleshchenkoFormFactors( -t1_, w1_, mx2 );
	- if ( inel_p2 ) fp2 = SzczurekUleshchenkoFormFactors( -t2_, w2_, my2 );
	+ fp1 = FormFactors::ProtonInelastic( cuts_.structure_functions, -t1_, w1_, mx2 );
	+ fp2 = FormFactors::ProtonInelastic( cuts_.structure_functions, -t2_, w2_, my2 );
	} break;
	}
	}

	inline bool
	GenericProcess::isKinematicsDefined()
	{
	// check the incoming state
	if ( !particles( Particle::IncomingBeam1 ).empty() && !particles( Particle::IncomingBeam2 ).empty() ) {
	is_incoming_state_set_ = true;
	}
	// check the outgoing state
	if ( !particles( Particle::OutgoingBeam1 ).empty()
	&& !particles( Particle::OutgoingBeam2 ).empty()
	&& !particles( Particle::CentralSystem ).empty() ) {
	is_outgoing_state_set_ = true;
	}
	// combine both states
	is_kinematics_set_ = is_incoming_state_set_ && is_outgoing_state_set_;
	return is_kinematics_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_physics_processes.cpp b/test/test_physics_processes.cpp
	index 74e1f90..a7cbad8 100644
	--- a/test/test_physics_processes.cpp
	+++ b/test/test_physics_processes.cpp
	@@ -1,123 +1,123 @@
	#include "CepGen/Generator.h"

	#include "CepGen/Processes/GamGamLL.h"
	#include "CepGen/Processes/PPtoLL.h"

	#include <assert.h>

	int
	main( int argc, char* argv[] )
	{
	typedef std::map<std::string,std::pair<double,double> > KinematicsMap;
	typedef std::map<float, KinematicsMap> ValuesAtCutMap;

	// values defined at pt(single lepton)>15 GeV, \|eta(single lepton)\|<2.5, mX<1000 GeV
	// process -> { pt cut -> { kinematics -> ( sigma, delta(sigma) ) } }
	std::map<std::string,ValuesAtCutMap> values_map = {
	//--- LPAIR values at sqrt(s) = 13 TeV
	{ "1_lpair", {
	{ 3.0, { // pt cut
	{ "1_elastic", { 2.0871703e1, 3.542e-2 } },
	{ "2_singlediss", { 1.5042536e1, 3.256e-2 } },
	{ "3_doublediss", { 1.38835e1, 4.03018e-2 } }
	} },
	{ 15.0, { // pt cut
	{ "1_elastic", { 4.1994803e-1, 8.328e-4 } },
	{ "2_singlediss", { 4.8504819e-1, 1.171e-3 } },
	{ "3_doublediss", { 6.35650e-1, 1.93968e-3 } }
	} },
	} },
	//--- PPTOLL values
	{ "2_pptoll", {
	{ 3.0, { // pt cut
	{ "1_elastic", { 2.0936541e1, 1.4096e-2 } },
	{ "2_singlediss", { 1.4844881e1, 2.0723e-2 } }, // SU, qt<50
	{ "3_doublediss", { 1.3580637e1, 2.2497e-2 } }, // SU, qt<50
	} },
	{ 15.0, { // pt cut
	{ "1_elastic", { 4.2515888e-1, 3.0351e-4 } },
	{ "2_singlediss", { 4.4903253e-1, 5.8970e-4 } }, // SU, qt<50
	{ "3_doublediss", { 5.1923819e-1, 9.6549e-4 } }, // SU, qt<50
	/*{ "2_singlediss", { 4.6710287e-1, 6.4842e-4 } }, // SU, qt<500
	{ "3_doublediss", { 5.6316353e-1, 1.1829e-3 } }, // SU, qt<500*/
	} },
	} },
	};

	const double num_sigma = 3.0;

	if ( argc == 1 \|\| strcmp( argv[1], "debug" ) != 0 ) {
	CepGen::Logger::get().level = CepGen::Logger::Nothing;
	}

	Timer tmr;
	CepGen::Generator mg;

	mg.parameters->kinematics.setSqrtS( 13.e3 );
	mg.parameters->kinematics.central_cuts[CepGen::Cuts::eta_single].in( -2.5, 2.5 );
	mg.parameters->kinematics.remnant_cuts[CepGen::Cuts::mass].max() = 1000.;
	//mg.parameters->vegas.ncvg = 50000;
	mg.parameters->vegas.itvg = 5;

	Information( Form( "Initial configuration time: %.3f ms", tmr.elapsed()*1.e3 ) );
	tmr.reset();

	unsigned short num_tests = 0, num_tests_passed = 0;
	std::vector<std::string> failed_tests;

	for ( const auto& values_vs_generator : values_map ) { // loop over all generators
	if ( values_vs_generator.first == "1_lpair" ) mg.parameters->setProcess( new CepGen::Process::GamGamLL );
	else if ( values_vs_generator.first == "2_pptoll" ) {
	mg.parameters->setProcess( new CepGen::Process::PPtoLL );
	- mg.parameters->kinematics.structure_functions = CepGen::SzczurekUleshchenko; //FIXME move to a dedicated class
	+ mg.parameters->kinematics.structure_functions = CepGen::StructureFunctionsType::SzczurekUleshchenko; //FIXME move to a dedicated class
	mg.parameters->kinematics.initial_cuts[CepGen::Cuts::qt].max() = 50.0;
	}
	else { InError( Form( "Unrecognized generator mode: %s", values_vs_generator.first.c_str() ) ); break; }

	for ( const auto& values_vs_cut : values_vs_generator.second ) { // loop over the single lepton pT cut
	mg.parameters->kinematics.central_cuts[CepGen::Cuts::pt_single].min() = values_vs_cut.first;
	for ( const auto& values_vs_kin : values_vs_cut.second ) { // loop over all possible kinematics
	if ( values_vs_kin.first == "1_elastic" ) mg.parameters->kinematics.mode = CepGen::Kinematics::ElasticElastic;
	else if ( values_vs_kin.first == "2_singlediss" ) mg.parameters->kinematics.mode = CepGen::Kinematics::InelasticElastic;
	else if ( values_vs_kin.first == "3_doublediss" ) mg.parameters->kinematics.mode = CepGen::Kinematics::InelasticInelastic;
	else { InError( Form( "Unrecognized kinematics mode: %s", values_vs_kin.first.c_str() ) ); break; }

	Information( Form( "Process: %s/%s\n\tConfiguration time: %.3f ms", values_vs_generator.first.c_str(), values_vs_kin.first.c_str(), tmr.elapsed()*1.e3 ) );
	tmr.reset();

	mg.clearRun();
	const double xsec_ref = values_vs_kin.second.first, err_xsec_ref = values_vs_kin.second.second;
	double xsec_cepgen, err_xsec_cepgen;
	mg.computeXsection( xsec_cepgen, err_xsec_cepgen );

	const double sigma = ( fabs( xsec_ref-xsec_cepgen ) ) / sqrt( err_xsec_cepgenerr_xsec_cepgen + err_xsec_referr_xsec_ref );

	Information( Form( "Computed cross section:\n\tRef. = %.3e +/- %.3e\n\tCepGen = %.3e +/- %.3e\n\tPull: %.6f", xsec_ref, err_xsec_ref, xsec_cepgen, err_xsec_cepgen, sigma ) );

	Information( Form( "Computation time: %.3f ms", tmr.elapsed()*1.e3 ) );
	tmr.reset();

	if ( fabs( sigma )<num_sigma ) num_tests_passed++;
	else {
	failed_tests.emplace_back( values_vs_generator.first+":"+
	Form( "pt-gt-%.1f", values_vs_cut.first )+":"+
	values_vs_kin.first+":"
	"ref="+Form( "%.3e", xsec_ref )+":"
	"got="+Form( "%.3e", xsec_cepgen )+":"
	"pull="+Form( "%.3f", sigma ) );
	}
	num_tests++;
	std::cout << "Test " << num_tests_passed << "/" << num_tests << " passed!" << std::endl;
	}
	}
	}
	if ( failed_tests.size() != 0 ) {
	std::ostringstream os;
	for ( const auto& fail : failed_tests ) os << " (*) " << fail << std::endl;
	throw CepGen::Exception( __PRETTY_FUNCTION__, Form( "Some tests failed!\n%s", os.str().c_str() ), CepGen::FatalError );
	}

	Information( "ALL TESTS PASSED!" );

	return 0;
	}
	diff --git a/test/utils/form_factors.cxx b/test/utils/form_factors.cxx
	index 4bd28fe..d332761 100644
	--- a/test/utils/form_factors.cxx
	+++ b/test/utils/form_factors.cxx
	@@ -1,90 +1,90 @@
	#include "CepGen/Physics/FormFactors.h"
	#include "CepGen/Physics/Particle.h"
	#include "test/Canvas.h"

	#include "TGraph.h"
	#include "TMultiGraph.h"

	#include <iostream>

	using namespace std;

	int
	main( int argc, char* argv[] )
	{
	const float min_q2 = 0.5, max_q2 = 1.e5, mx2 = ( argc>1 ) ? atof( argv[1] ) : 1.e4;
	const char* mx2_str = ( argc>2 ) ? argv[2] : "10^{4}";
	const unsigned int npoints = 1000;

	TGraph g_sy_fe_100, g_sy_fm_100;
	//TGraph g_fb_fe_100, g_fb_fm_100;
	TGraph g_su_fe_100, g_su_fm_100;

	const float mp2 = pow( CepGen::Particle::massFromPDGId( CepGen::Particle::Proton ), 2 );

	for ( unsigned int i=0; i<npoints; i++ ) {
	const float q2 = min_q2 + i*( max_q2-min_q2 )/(npoints-1);

	- CepGen::FormFactors ff_sy = CepGen::SuriYennieFormFactors( q2, mp2, mx2 );
	+ CepGen::FormFactors ff_sy = CepGen::FormFactors::SuriYennie( q2, mp2, mx2 );
	g_sy_fe_100.SetPoint( i, q2, ff_sy.FE );
	g_sy_fm_100.SetPoint( i, q2, ff_sy.FM );

	/*CepGen::FormFactors ff_fb = CepGen::FioreBrasseFormFactors( q2, mp2, mx2 );
	g_fb_fe_100.SetPoint( i, q2, ff_fb.FE );
	g_fb_fm_100.SetPoint( i, q2, ff_fb.FM );*/

	- CepGen::FormFactors ff_su = CepGen::SzczurekUleshchenkoFormFactors( q2, mp2, mx2 );
	+ CepGen::FormFactors ff_su = CepGen::FormFactors::SzczurekUleshchenko( q2, mp2, mx2 );
	g_su_fe_100.SetPoint( i, q2, ff_su.FE );
	g_su_fm_100.SetPoint( i, q2, ff_su.FM );
	cout << q2 << "\t" << ff_su.FM << endl;
	}

	CepGen::Canvas c( "test", Form( "CepGen proton form factors, M_{X}^{2} = %s GeV^{2}", mx2_str ) );
	c.SetLegendX1( 0.4 );

	TMultiGraph mg;

	g_sy_fe_100.SetLineWidth( 3 );
	mg.Add( &g_sy_fe_100, "l" );
	c.AddLegendEntry( &g_sy_fe_100, "Suri-Yennie, F_{E}", "l" );

	g_sy_fm_100.SetLineStyle( 2 );
	g_sy_fm_100.SetLineWidth( 3 );
	mg.Add( &g_sy_fm_100, "l" );
	c.AddLegendEntry( &g_sy_fm_100, "Suri-Yennie, F_{M}", "l" );

	/*g_fb_fe_100.SetLineColor( kRed+1 );
	g_fb_fe_100.SetLineWidth( 3 );
	mg.Add( &g_fb_fe_100, "l" );
	c.AddLegendEntry( &g_fb_fe_100, "Fiore-Brasse, F_{E}", "l" );

	g_fb_fm_100.SetLineStyle( 2 );
	g_fb_fm_100.SetLineColor( kRed+1 );
	g_fb_fm_100.SetLineWidth( 3 );
	mg.Add( &g_fb_fm_100, "l" );
	c.AddLegendEntry( &g_fb_fm_100, "Fiore-Brasse, F_{M}", "l" );*/

	g_su_fe_100.SetLineColor( kGreen+2 );
	g_su_fe_100.SetLineWidth( 3 );
	mg.Add( &g_su_fe_100, "l" );
	c.AddLegendEntry( &g_su_fe_100, "Szczurek-Uleshchenko, F_{E}", "l" );

	g_su_fm_100.SetLineStyle( 2 );
	g_su_fm_100.SetLineColor( kGreen+2 );
	g_su_fm_100.SetLineWidth( 3 );
	mg.Add( &g_su_fm_100, "l" );
	c.AddLegendEntry( &g_su_fm_100, "Szczurek-Uleshchenko, F_{M}", "l" );

	mg.Draw( "alpr" );
	mg.SetTitle( "Q^{2} (GeV^{2})\\Proton form factor" );

	c.Prettify( mg.GetHistogram() );
	mg.GetYaxis()->SetRangeUser( 1.e-10, 10. );
	mg.GetXaxis()->SetLimits( min_q2, max_q2 );
	c.SetLogx();
	c.SetLogy();

	c.Save( "pdf" );

	return 0;
	}
	diff --git a/test/utils/form_factors.cxx b/test/utils/form_factors_vs_xbj.cxx
	similarity index 56%
	copy from test/utils/form_factors.cxx
	copy to test/utils/form_factors_vs_xbj.cxx
	index 4bd28fe..67fe939 100644
	--- a/test/utils/form_factors.cxx
	+++ b/test/utils/form_factors_vs_xbj.cxx
	@@ -1,90 +1,91 @@
	#include "CepGen/Physics/FormFactors.h"
	#include "CepGen/Physics/Particle.h"
	#include "test/Canvas.h"

	#include "TGraph.h"
	#include "TMultiGraph.h"

	#include <iostream>

	using namespace std;

	int
	main( int argc, char* argv[] )
	{
	- const float min_q2 = 0.5, max_q2 = 1.e5, mx2 = ( argc>1 ) ? atof( argv[1] ) : 1.e4;
	- const char* mx2_str = ( argc>2 ) ? argv[2] : "10^{4}";
	+ const float min_xbj = 1.e-3, max_xbj = 1.0, q2 = ( argc>1 ) ? atof( argv[1] ) : 2.5;
	+ const char* q2_str = ( argc>2 ) ? argv[2] : "2.5";
	const unsigned int npoints = 1000;

	TGraph g_sy_fe_100, g_sy_fm_100;
	- //TGraph g_fb_fe_100, g_fb_fm_100;
	+ TGraph g_fb_fe_100, g_fb_fm_100;
	TGraph g_su_fe_100, g_su_fm_100;

	const float mp2 = pow( CepGen::Particle::massFromPDGId( CepGen::Particle::Proton ), 2 );

	for ( unsigned int i=0; i<npoints; i++ ) {
	- const float q2 = min_q2 + i*( max_q2-min_q2 )/(npoints-1);
	+ const float xbj = min_xbj + i*( max_xbj-min_xbj )/( npoints-1 );
	+ const float mx2 = mp2 + ( 1.-xbj )*q2;

	- CepGen::FormFactors ff_sy = CepGen::SuriYennieFormFactors( q2, mp2, mx2 );
	- g_sy_fe_100.SetPoint( i, q2, ff_sy.FE );
	- g_sy_fm_100.SetPoint( i, q2, ff_sy.FM );
	+ CepGen::FormFactors ff_sy = CepGen::FormFactors::SuriYennie( q2, mp2, mx2 );
	+ g_sy_fe_100.SetPoint( i, xbj, ff_sy.FE );
	+ g_sy_fm_100.SetPoint( i, xbj, ff_sy.FM );

	- /*CepGen::FormFactors ff_fb = CepGen::FioreBrasseFormFactors( q2, mp2, mx2 );
	- g_fb_fe_100.SetPoint( i, q2, ff_fb.FE );
	- g_fb_fm_100.SetPoint( i, q2, ff_fb.FM );*/
	+ CepGen::FormFactors ff_fb = CepGen::FormFactors::FioreBrasse( q2, mp2, mx2 );
	+ g_fb_fe_100.SetPoint( i, xbj, ff_fb.FE );
	+ g_fb_fm_100.SetPoint( i, xbj, ff_fb.FM );

	- CepGen::FormFactors ff_su = CepGen::SzczurekUleshchenkoFormFactors( q2, mp2, mx2 );
	- g_su_fe_100.SetPoint( i, q2, ff_su.FE );
	- g_su_fm_100.SetPoint( i, q2, ff_su.FM );
	- cout << q2 << "\t" << ff_su.FM << endl;
	+ CepGen::FormFactors ff_su = CepGen::FormFactors::SzczurekUleshchenko( q2, mp2, mx2 );
	+ g_su_fe_100.SetPoint( i, xbj, ff_su.FE );
	+ g_su_fm_100.SetPoint( i, xbj, ff_su.FM );
	}

	- CepGen::Canvas c( "test", Form( "CepGen proton form factors, M_{X}^{2} = %s GeV^{2}", mx2_str ) );
	+ CepGen::Canvas c( "test", Form( "CepGen proton form factors, Q^{2} = %s GeV^{2}", q2_str ) );
	c.SetLegendX1( 0.4 );

	TMultiGraph mg;

	g_sy_fe_100.SetLineWidth( 3 );
	mg.Add( &g_sy_fe_100, "l" );
	c.AddLegendEntry( &g_sy_fe_100, "Suri-Yennie, F_{E}", "l" );

	g_sy_fm_100.SetLineStyle( 2 );
	g_sy_fm_100.SetLineWidth( 3 );
	mg.Add( &g_sy_fm_100, "l" );
	c.AddLegendEntry( &g_sy_fm_100, "Suri-Yennie, F_{M}", "l" );

	- /*g_fb_fe_100.SetLineColor( kRed+1 );
	+ g_fb_fe_100.SetLineColor( kRed+1 );
	g_fb_fe_100.SetLineWidth( 3 );
	mg.Add( &g_fb_fe_100, "l" );
	c.AddLegendEntry( &g_fb_fe_100, "Fiore-Brasse, F_{E}", "l" );

	g_fb_fm_100.SetLineStyle( 2 );
	g_fb_fm_100.SetLineColor( kRed+1 );
	g_fb_fm_100.SetLineWidth( 3 );
	mg.Add( &g_fb_fm_100, "l" );
	- c.AddLegendEntry( &g_fb_fm_100, "Fiore-Brasse, F_{M}", "l" );*/
	+ c.AddLegendEntry( &g_fb_fm_100, "Fiore-Brasse, F_{M}", "l" );

	g_su_fe_100.SetLineColor( kGreen+2 );
	g_su_fe_100.SetLineWidth( 3 );
	mg.Add( &g_su_fe_100, "l" );
	c.AddLegendEntry( &g_su_fe_100, "Szczurek-Uleshchenko, F_{E}", "l" );

	g_su_fm_100.SetLineStyle( 2 );
	g_su_fm_100.SetLineColor( kGreen+2 );
	g_su_fm_100.SetLineWidth( 3 );
	mg.Add( &g_su_fm_100, "l" );
	c.AddLegendEntry( &g_su_fm_100, "Szczurek-Uleshchenko, F_{M}", "l" );

	mg.Draw( "alpr" );
	- mg.SetTitle( "Q^{2} (GeV^{2})\\Proton form factor" );
	+ mg.SetTitle( "x_{Bj}\\Proton form factor" );

	c.Prettify( mg.GetHistogram() );
	- mg.GetYaxis()->SetRangeUser( 1.e-10, 10. );
	- mg.GetXaxis()->SetLimits( min_q2, max_q2 );
	- c.SetLogx();
	- c.SetLogy();
	+ //mg.GetYaxis()->SetRangeUser( 1.e-10, 10. );
	+ //mg.GetYaxis()->SetRangeUser( 0., 0.6 );
	+ mg.GetXaxis()->SetLimits( min_xbj, max_xbj );
	+ //c.SetLogx();
	+ //c.SetLogy();

	c.Save( "pdf" );

	return 0;
	}
	diff --git a/test/utils/structure_functions.cxx b/test/utils/structure_functions.cxx
	new file mode 100644
	index 0000000..5a16ed2
	--- /dev/null
	+++ b/test/utils/structure_functions.cxx
	@@ -0,0 +1,97 @@
	+#include "CepGen/Physics/FormFactors.h"
	+#include "CepGen/Physics/Particle.h"
	+#include "test/Canvas.h"
	+
	+#include "TGraph.h"
	+#include "TMultiGraph.h"
	+
	+#include <iostream>
	+
	+using namespace std;
	+
	+int
	+main( int argc, char* argv[] )
	+{
	+ const float min_xbj = 1.e-3, max_xbj = 0.99, q2 = ( argc>1 ) ? atof( argv[1] ) : 2.5;
	+ const char* q2_str = ( argc>2 ) ? argv[2] : std::to_string( q2 ).c_str();
	+ const unsigned int npoints = 5000;
	+
	+ TGraph g_sy_f1, g_sy_f2;
	+ TGraph g_fb_f1, g_fb_f2;
	+ TGraph g_su_f1, g_su_f2;
	+
	+ const bool use_logarithmic_x = false;
	+
	+ for ( unsigned int i=0; i<npoints; i++ ) {
	+ float xbj;
	+ if ( use_logarithmic_x ) {
	+ const float min_lxbj = log10( min_xbj ), max_lxbj = log10( max_xbj );
	+ xbj = pow( 10, min_lxbj + i*( max_lxbj-min_lxbj )/( npoints-1 ) );
	+ std::cout << min_lxbj << "\t" << max_lxbj << "\t" << xbj << std::endl;
	+ }
	+ else xbj = min_xbj + i*( max_xbj-min_xbj )/( npoints-1 );
	+
	+ /*CepGen::StructureFunctions sf_sy = CepGen::StructureFunctions::SuriYennie( q2, mp2, mx2 );
	+ g_sy_f1.SetPoint( i, xbj, ff_sy.F1 );
	+ g_sy_f2.SetPoint( i, xbj, ff_sy.F2 );*/
	+
	+ CepGen::StructureFunctions sf_fb = CepGen::StructureFunctions::FioreBrasse( q2, xbj );
	+ g_fb_f1.SetPoint( i, xbj, sf_fb.F1 );
	+ g_fb_f2.SetPoint( i, xbj, sf_fb.F2*xbj );
	+
	+ CepGen::StructureFunctions sf_su = CepGen::StructureFunctions::SzczurekUleshchenko( q2, xbj );
	+ g_su_f1.SetPoint( i, xbj, sf_su.F1 );
	+ g_su_f2.SetPoint( i, xbj, sf_su.F2*xbj );
	+ std::cout << sf_fb << std::endl;
	+ }
	+
	+ CepGen::Canvas c( "test", Form( "CepGen proton structure functions, Q^{2} = %s GeV^{2}", q2_str ) );
	+ c.SetLegendX1( 0.4 );
	+
	+ TMultiGraph mg;
	+
	+ /*g_sy_f1.SetLineWidth( 3 );
	+ mg.Add( &g_sy_f1, "l" );
	+ c.AddLegendEntry( &g_sy_f1, "Suri-Yennie, F_{E}", "l" );
	+
	+ g_sy_f2.SetLineStyle( 2 );
	+ g_sy_f2.SetLineWidth( 3 );
	+ mg.Add( &g_sy_f2, "l" );
	+ c.AddLegendEntry( &g_sy_f2, "Suri-Yennie, F_{M}", "l" );*/
	+
	+ /*g_fb_f1.SetLineColor( kRed+1 );
	+ g_fb_f1.SetLineWidth( 3 );
	+ mg.Add( &g_fb_f1, "l" );
	+ c.AddLegendEntry( &g_fb_f1, "Fiore-Brasse, F_{1}", "l" );*/
	+
	+ //g_fb_f2.SetLineStyle( 2 );
	+ g_fb_f2.SetLineColor( kRed+1 );
	+ g_fb_f2.SetLineWidth( 3 );
	+ mg.Add( &g_fb_f2, "l" );
	+ c.AddLegendEntry( &g_fb_f2, "Fiore-Brasse, F_{2}", "l" );
	+
	+ /*g_su_f1.SetLineColor( kGreen+2 );
	+ g_su_f1.SetLineWidth( 3 );
	+ mg.Add( &g_su_f1, "l" );
	+ c.AddLegendEntry( &g_su_f1, "Szczurek-Uleshchenko, F_{1}", "l" );*/
	+
	+ //g_su_f2.SetLineStyle( 2 );
	+ g_su_f2.SetLineColor( kGreen+2 );
	+ g_su_f2.SetLineWidth( 3 );
	+ mg.Add( &g_su_f2, "l" );
	+ c.AddLegendEntry( &g_su_f2, "Szczurek-Uleshchenko, F_{2}", "l" );
	+
	+ mg.Draw( "alpr" );
	+ mg.SetTitle( "x_{Bj}\\Proton form factor" );
	+
	+ c.Prettify( mg.GetHistogram() );
	+ //mg.GetYaxis()->SetRangeUser( 1.e-10, 10. );
	+ //mg.GetYaxis()->SetRangeUser( 0., 0.6 );
	+ mg.GetXaxis()->SetLimits( min_xbj, max_xbj );
	+ if ( use_logarithmic_x ) c.SetLogx();
	+ //c.SetLogy();
	+
	+ c.Save( "pdf" );
	+
	+ return 0;
	+}

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 8:35 PM (1 d, 3 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3806066
Default Alt Text: (73 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions