diff --git a/CepGen/Cards/ConfigHandler.cpp b/CepGen/Cards/ConfigHandler.cpp
index 1dcf4f2..cbff4fe 100644
--- a/CepGen/Cards/ConfigHandler.cpp
+++ b/CepGen/Cards/ConfigHandler.cpp
@@ -1,257 +1,257 @@
#include "ConfigHandler.h"
#ifdef LIBCONFIG
namespace CepGen
{
namespace Cards
{
//----- specialization for CepGen input cards
ConfigHandler::ConfigHandler( const char* file )
{
#ifdef LIBCONFIG
libconfig::Config cfg;
try { cfg.readFile( file ); } catch ( const libconfig::ParseException& pe ) {
FatalError( Form( "Failed to parse the configuration card \"%s\".\n\tParser error: %s (at line %d)", file, pe.getError(), pe.getLine() ) );
}
try {
const libconfig::Setting& root = cfg.getRoot();
const libconfig::Setting& proc = root["process"];
//--- type of process to consider
const char* proc_name = proc["name"]; const std::string str_proc_name = proc_name;
if ( str_proc_name == "lpair" ) params_.setProcess( new Process::GamGamLL );
else if ( str_proc_name == "pptoll" ) params_.setProcess( new Process::PPtoLL );
else if ( str_proc_name == "pptoww" ) params_.setProcess( new Process::PPtoWW );
else FatalError( Form( "Unrecognised process: %s", proc_name ) );
//--- process mode
int int_mode; const char* mode;
if ( proc.lookupValue( "mode", int_mode ) ) {
params_.kinematics.mode = (Kinematics::ProcessMode)int_mode;
}
else if ( proc.lookupValue( "mode", mode ) ) {
const std::string str_mode = mode;
if ( str_mode == "elastic/elastic" ) params_.kinematics.mode = Kinematics::ProcessMode::ElasticElastic;
else if ( str_mode == "elastic/inelastic" ) params_.kinematics.mode = Kinematics::ProcessMode::ElasticInelastic;
else if ( str_mode == "inelastic/elastic" ) params_.kinematics.mode = Kinematics::ProcessMode::InelasticElastic;
else if ( str_mode == "inelastic/inelastic" ) params_.kinematics.mode = Kinematics::ProcessMode::InelasticInelastic;
else FatalError( Form( "Unrecognised interaction mode: %s", mode ) );
}
//--- process kinematics
if ( proc.exists( "in_kinematics" ) ) parseIncomingKinematics( proc["in_kinematics"] );
if ( proc.exists( "out_kinematics" ) ) parseOutgoingKinematics( proc["out_kinematics"] );
//--- generation parameters
if ( root.exists( "vegas" ) ) parseVegas( root["vegas"] );
if ( root.exists( "generator" ) ) parseGenerator( root["generator"] );
//--- taming functions
if ( proc.exists( "taming_functions" ) ) parseTamingFunctions( proc["taming_functions"] );
} catch ( const libconfig::SettingNotFoundException& nfe ) {
FatalError( Form( "Failed to retrieve the field \"%s\".", nfe.getPath() ) );
} catch ( const libconfig::SettingTypeException& te ) {
FatalError( Form( "Field \"%s\" has wrong type.", te.getPath() ) );
}
#else
InWarning( "libconfig++ is not present on this machine" );
#endif
}
#ifdef LIBCONFIG
void
ConfigHandler::parseIncomingKinematics( const libconfig::Setting& kin )
{
try {
if ( kin.exists( "beam1_pz" ) ) params_.kinematics.inp.first = (double)kin["beam1_pz"];
if ( kin.exists( "beam2_pz" ) ) params_.kinematics.inp.second = (double)kin["beam2_pz"];
if ( kin.exists( "structure_functions" ) ) {
const char* sf = kin["structure_functions" ]; const std::string sf_str = sf;
if ( sf_str == "electron" ) params_.kinematics.structure_functions = Electron;
else if ( sf_str == "elastic proton" ) params_.kinematics.structure_functions = ElasticProton;
else if ( sf_str == "Suri-Yennie" ) params_.kinematics.structure_functions = SuriYennie;
else if ( sf_str == "Suri-Yennie;lowQ2" ) params_.kinematics.structure_functions = SuriYennieLowQ2;
else if ( sf_str == "Szczurek-Uleshchenko" ) params_.kinematics.structure_functions = SzczurekUleshchenko;
- else if ( sf_str == "Fiore;valence" ) params_.kinematics.structure_functions = FioreVal;
- else if ( sf_str == "Fiore;sea" ) params_.kinematics.structure_functions = FioreSea;
- else if ( sf_str == "Fiore" ) params_.kinematics.structure_functions = Fiore;
+ else if ( sf_str == "Fiore" ) params_.kinematics.structure_functions = FioreBrasse;
else if ( sf_str == "ALLM" ) params_.kinematics.structure_functions = ALLM97;
else if ( sf_str == "ALLM;91" ) params_.kinematics.structure_functions = ALLM91;
else if ( sf_str == "ALLM;97" ) params_.kinematics.structure_functions = ALLM97;
else if ( sf_str == "ALLM;HHT" ) params_.kinematics.structure_functions = ALLM_HHT;
else if ( sf_str == "ALLM;HHT-FT" ) params_.kinematics.structure_functions = ALLM_HHT_FT;
+ else if ( sf_str == "ALLM;GD07p" ) params_.kinematics.structure_functions = ALLM_GD07p;
+ else if ( sf_str == "ALLM;GD11p" ) params_.kinematics.structure_functions = ALLM_GD11p;
else FatalError( Form( "Invalid structure functions mode: %s", sf ) );
}
} catch ( const libconfig::SettingNotFoundException& nfe ) {
FatalError( Form( "Failed to retrieve the field \"%s\".", nfe.getPath() ) );
} catch ( const libconfig::SettingTypeException& te ) {
FatalError( Form( "Field \"%s\" has wrong type.", te.getPath() ) );
}
}
void
ConfigHandler::parseOutgoingKinematics( const libconfig::Setting& kin )
{
try {
if ( kin.exists( "pair" ) ) {
Particle::ParticleCode pair = (Particle::ParticleCode)(int)kin["pair"];
params_.kinematics.central_system = { pair, pair };
}
if ( kin.exists( "min_pt" ) ) params_.kinematics.central_cuts[Cuts::pt_single].min() = (double)kin["min_pt"];
if ( kin.exists( "max_pt" ) ) params_.kinematics.central_cuts[Cuts::pt_single].max() = (double)kin["max_pt"];
if ( kin.exists( "min_ptdiff" ) ) params_.kinematics.central_cuts[Cuts::pt_diff].min() = (double)kin["min_ptdiff"];
if ( kin.exists( "max_ptdiff" ) ) params_.kinematics.central_cuts[Cuts::pt_diff].max() = (double)kin["max_ptdiff"];
if ( kin.exists( "min_dely" ) ) params_.kinematics.central_cuts[Cuts::dely].min() = (double)kin["min_dely"];
if ( kin.exists( "max_dely" ) ) params_.kinematics.central_cuts[Cuts::dely].max() = (double)kin["max_dely"];
if ( kin.exists( "min_energy" ) ) params_.kinematics.central_cuts[Cuts::energy_single].min() = (double)kin["min_energy"];
if ( kin.exists( "max_energy" ) ) params_.kinematics.central_cuts[Cuts::energy_single].max() = (double)kin["max_energy"];
if ( kin.exists( "min_eta" ) ) params_.kinematics.central_cuts[Cuts::eta_single].min() = (double)kin["min_eta"];
if ( kin.exists( "max_eta" ) ) params_.kinematics.central_cuts[Cuts::eta_single].max() = (double)kin["max_eta"];
if ( kin.exists( "min_rapidity" ) ) params_.kinematics.central_cuts[Cuts::rapidity_single].min() = (double)kin["min_rapidity"];
if ( kin.exists( "max_rapidity" ) ) params_.kinematics.central_cuts[Cuts::rapidity_single].max() = (double)kin["max_rapidity"];
if ( kin.exists( "min_mx" ) ) params_.kinematics.remnant_cuts[Cuts::mass].min() = (double)kin["min_mx"];
if ( kin.exists( "max_mx" ) ) params_.kinematics.remnant_cuts[Cuts::mass].max() = (double)kin["max_mx"];
} catch ( const libconfig::SettingNotFoundException& nfe ) {
FatalError( Form( "Failed to retrieve the field \"%s\".", nfe.getPath() ) );
} catch ( const libconfig::SettingTypeException& te ) {
FatalError( Form( "Field \"%s\" has wrong type.", te.getPath() ) );
}
}
void
ConfigHandler::parseVegas( const libconfig::Setting& veg )
{
try {
if ( veg.exists( "num_points" ) ) params_.vegas.npoints = (int)veg["num_points"];
if ( veg.exists( "num_integration_calls" ) ) params_.vegas.ncvg = (int)veg["num_integration_calls"];
if ( veg.exists( "num_integration_iterations" ) ) params_.vegas.itvg = (int)veg["num_integration_iterations"];
if ( veg.exists( "seed" ) ) params_.vegas.seed = (unsigned long)veg["seed"];
} catch ( const libconfig::SettingNotFoundException& nfe ) {
FatalError( Form( "Failed to retrieve the field \"%s\".", nfe.getPath() ) );
}
}
void
ConfigHandler::parseGenerator( const libconfig::Setting& gen )
{
params_.generation.enabled = true;
try {
if ( gen.exists( "num_events" ) ) params_.generation.maxgen = (int)gen["num_events"];
if ( gen.exists( "print_every" ) ) params_.generation.gen_print_every = (int)gen["print_every"];
} catch ( const libconfig::SettingNotFoundException& nfe ) {
FatalError( Form( "Failed to retrieve the field \"%s\".", nfe.getPath() ) );
}
}
void
ConfigHandler::parseTamingFunctions( const libconfig::Setting& tf )
{
if ( !tf.isList() ) FatalError( "The taming functions definition must be wrapped within a list!" );
for ( unsigned short i = 0; i < tf.getLength(); ++i ) {
params_.taming_functions.add( tf[i]["variable"], tf[i]["expression"] );
}
}
void
ConfigHandler::writeProcess( const Parameters* params, libconfig::Setting& root )
{
libconfig::Setting& proc = root.add( "process", libconfig::Setting::TypeGroup );
proc.add( "name", libconfig::Setting::TypeString ) = params->processName();
std::ostringstream os; os << params->kinematics.mode;
proc.add( "mode", libconfig::Setting::TypeString ) = os.str();
}
void
ConfigHandler::writeIncomingKinematics( const Parameters* params, libconfig::Setting& root )
{
libconfig::Setting& kin = root.add( "in_kinematics", libconfig::Setting::TypeGroup );
kin.add( "beam1_pz", libconfig::Setting::TypeFloat ) = params->kinematics.inp.first;
kin.add( "beam2_pz", libconfig::Setting::TypeFloat ) = params->kinematics.inp.second;
std::ostringstream os; os << params->kinematics.structure_functions;
kin.add( "structure_function", libconfig::Setting::TypeString ) = os.str();
}
void
ConfigHandler::writeOutgoingKinematics( const Parameters* params, libconfig::Setting& root )
{
libconfig::Setting& kin = root.add( "out_kinematics", libconfig::Setting::TypeGroup );
if ( params->kinematics.central_system.size() > 0 )
kin.add( "pair", libconfig::Setting::TypeInt ) = (int)params->kinematics.central_system[0];
if ( params->kinematics.central_cuts.count( Cuts::pt_single ) ) {
kin.add( "min_pt", libconfig::Setting::TypeFloat ) = params->kinematics.central_cuts.at( Cuts::pt_single ).min();
kin.add( "max_pt", libconfig::Setting::TypeFloat ) = params->kinematics.central_cuts.at( Cuts::pt_single ).max();
}
if ( params->kinematics.central_cuts.count( Cuts::pt_diff ) ) {
kin.add( "min_ptdiff", libconfig::Setting::TypeFloat ) = params->kinematics.central_cuts.at( Cuts::pt_diff ).min();
kin.add( "max_ptdiff", libconfig::Setting::TypeFloat ) = params->kinematics.central_cuts.at( Cuts::pt_diff ).max();
}
if ( params->kinematics.central_cuts.count( Cuts::dely ) ) {
kin.add( "min_dely", libconfig::Setting::TypeFloat ) = params->kinematics.central_cuts.at( Cuts::dely ).min();
kin.add( "max_dely", libconfig::Setting::TypeFloat ) = params->kinematics.central_cuts.at( Cuts::dely ).max();
}
if ( params->kinematics.central_cuts.count( Cuts::energy_single ) ) {
kin.add( "min_energy", libconfig::Setting::TypeFloat ) = params->kinematics.central_cuts.at( Cuts::energy_single ).min();
kin.add( "max_energy", libconfig::Setting::TypeFloat ) = params->kinematics.central_cuts.at( Cuts::energy_single ).max();
}
if ( params->kinematics.central_cuts.count( Cuts::eta_single ) ) {
kin.add( "min_eta", libconfig::Setting::TypeFloat ) = params->kinematics.central_cuts.at( Cuts::eta_single ).min();
kin.add( "max_eta", libconfig::Setting::TypeFloat ) = params->kinematics.central_cuts.at( Cuts::eta_single ).max();
}
if ( params->kinematics.remnant_cuts.count( Cuts::mass ) ) {
kin.add( "min_mx", libconfig::Setting::TypeFloat ) = params->kinematics.remnant_cuts.at( Cuts::mass ).min();
kin.add( "max_mx", libconfig::Setting::TypeFloat ) = params->kinematics.remnant_cuts.at( Cuts::mass ).max();
}
}
void
ConfigHandler::writeTamingFunctions( const Parameters* params, libconfig::Setting& root )
{
libconfig::Setting& tf = root.add( "taming_functions", libconfig::Setting::TypeList );
for ( std::map<std::string,TamingFunction>::const_iterator it = params->taming_functions.begin(); it != params->taming_functions.end(); ++it ) {
libconfig::Setting& fun = tf.add( libconfig::Setting::TypeGroup );
fun.add( "variable", libconfig::Setting::TypeString ) = it->first;
fun.add( "expression", libconfig::Setting::TypeString ) = it->second.expression;
}
}
void
ConfigHandler::writeVegas( const Parameters* params, libconfig::Setting& root )
{
libconfig::Setting& veg = root.add( "vegas", libconfig::Setting::TypeGroup );
veg.add( "num_points", libconfig::Setting::TypeInt ) = (int)params->vegas.npoints;
veg.add( "num_integration_calls", libconfig::Setting::TypeInt ) = (int)params->vegas.ncvg;
veg.add( "num_integration_iterations", libconfig::Setting::TypeInt ) = (int)params->vegas.itvg;
veg.add( "seed", libconfig::Setting::TypeInt64 ) = (long)params->vegas.seed;
}
void
ConfigHandler::writeGenerator( const Parameters* params, libconfig::Setting& root )
{
if ( !params->generation.enabled ) return;
libconfig::Setting& gen = root.add( "generator", libconfig::Setting::TypeGroup );
gen.add( "num_events", libconfig::Setting::TypeInt ) = (int)params->generation.maxgen;
gen.add( "print_every", libconfig::Setting::TypeInt ) = (int)params->generation.gen_print_every;
}
#endif
void
ConfigHandler::store( const Parameters* params, const char* file )
{
#ifdef LIBCONFIG
libconfig::Config cfg;
libconfig::Setting& root = cfg.getRoot();
writeProcess( params, root );
writeIncomingKinematics( params, root["process"] );
writeOutgoingKinematics( params, root["process"] );
writeTamingFunctions( params, root["process"] );
writeVegas( params, root );
writeGenerator( params, root );
cfg.writeFile( file );
#endif
}
}
}
#endif
diff --git a/CepGen/Physics/FormFactors.cpp b/CepGen/Physics/FormFactors.cpp
index e27b879..905adb0 100644
--- a/CepGen/Physics/FormFactors.cpp
+++ b/CepGen/Physics/FormFactors.cpp
@@ -1,74 +1,72 @@
#include "FormFactors.h"
namespace CepGen
{
FormFactors
FormFactors::Trivial()
{
return FormFactors( 1.0, 1.0 );
}
FormFactors
FormFactors::ProtonElastic( double q2 )
{
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.79*GE, GM2 = GM*GM;
return FormFactors( ( 4.*mp2*GE2 + q2*GM2 ) / ( 4.*mp2 + q2 ), GM2 );
}
FormFactors
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:
+ case StructureFunctionsType::FioreBrasse:
return FormFactors::FioreBrasse( q2, mi2, mf2 );
default: throw Exception( __PRETTY_FUNCTION__, "Invalid structure functions required!", FatalError );
}
}
FormFactors
FormFactors::SuriYennie( double q2, double mi2, double mf2 )
{
const double x = q2 / ( q2 + mf2 - mi2 );
SF::SuriYennie suriyennie;
const StructureFunctions sy = suriyennie( q2, x );
//std::cout << "---> " << sy.FM << "\t" << sy.F2*x/q2 << "\t" << sy.F2*x*sqrt(mi2)/q2 << std::endl;
return FormFactors( sy.F2 * x * sqrt( mi2 ) / q2, sy.FM ); //FIXME
}
FormFactors
FormFactors::FioreBrasse( double q2, double mi2, double mf2 )
{
const double x = q2 / ( q2 + mf2 - mi2 );
SF::FioreBrasse fb;
StructureFunctions sf = fb( q2, x );
return FormFactors( sf.F2 * x / q2, -2.*sf.F1 / q2 );
}
FormFactors
FormFactors::SzczurekUleshchenko( double q2, double mi2, double mf2 )
{
const double x = q2 / ( q2 + mf2 - mi2 );
SF::SzczurekUleshchenko su;
StructureFunctions sf = su( q2, x );
return FormFactors( sf.F2 * x / q2, -2.*sf.F1 / q2 );
}
std::ostream&
operator<<( std::ostream& os, const FormFactors& ff )
{
os << Form( "Form factors: electric: Fe = %.3e ; magnetic: Fm = %.3e", ff.FE, ff.FM ).c_str();
return os;
}
}
diff --git a/CepGen/Processes/GenericKTProcess.cpp b/CepGen/Processes/GenericKTProcess.cpp
index dabab26..1bee34b 100644
--- a/CepGen/Processes/GenericKTProcess.cpp
+++ b/CepGen/Processes/GenericKTProcess.cpp
@@ -1,279 +1,285 @@
#include "GenericKTProcess.h"
namespace CepGen
{
namespace Process
{
GenericKTProcess::GenericKTProcess( const std::string& name,
const std::string& description,
const unsigned int& num_user_dimensions,
const std::array<Particle::ParticleCode,2>& partons,
const std::vector<Particle::ParticleCode>& central ) :
GenericProcess( name, description+" (kT-factorisation approach)" ),
kNumUserDimensions( num_user_dimensions ),
kIntermediateParts( partons ), kProducedParts( central )
{}
GenericKTProcess::~GenericKTProcess()
{}
void
GenericKTProcess::addEventContent()
{
GenericProcess::setEventContent(
{ // incoming state
{ Particle::IncomingBeam1, Particle::Proton },
{ Particle::IncomingBeam2, Particle::Proton },
{ Particle::Parton1, kIntermediateParts[0] },
{ Particle::Parton2, kIntermediateParts[1] }
},
{ // outgoing state
{ Particle::OutgoingBeam1, { Particle::Proton } },
{ Particle::OutgoingBeam2, { Particle::Proton } },
{ Particle::CentralSystem, kProducedParts }
}
);
setExtraContent();
}
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_PI*x( 2 );
phi_qt2_ = 2.*M_PI*x( 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_ = elasticFlux( x1, q1t2 );
flux2_ = elasticFlux( x2, q2t2 );
break;
case Kinematics::ElasticInelastic:
flux1_ = elasticFlux( x1, q1t2 );
flux2_ = inelasticFlux( x2, q2t2, MY_ );
break;
case Kinematics::InelasticElastic:
flux1_ = inelasticFlux( x1, q1t2, MX_ );
flux2_ = elasticFlux( x2, q2t2 );
break;
case Kinematics::InelasticInelastic:
flux1_ = inelasticFlux( x1, q1t2, MX_ );
flux2_ = inelasticFlux( x2, q2t2, MY_ );
break;
default: return;
}
flux1_ = std::max( flux1_, 1.e-20 );
flux2_ = std::max( flux2_, 1.e-20 );
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 );
op1.setMomentum( PX_ );
op2.setMomentum( PY_ );
switch ( cuts_.mode ) {
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( MY_ );
break;
case Kinematics::InelasticElastic:
op1.setStatus( Particle::Undecayed ); op1.setMass( MX_ );
op2.setStatus( Particle::FinalState );
break;
case Kinematics::InelasticInelastic:
op1.setStatus( Particle::Undecayed ); op1.setMass( MX_ );
op2.setStatus( Particle::Undecayed ); op2.setMass( MY_ );
break;
default: { FatalError( "This kT factorisation process is intended for p-on-p collisions! Aborting!" ); } break;
}
//=================================================================
// 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_min = x*x*mp2/( 1.-x ), Q2_ela = Q2_min + kt2/( 1.-x );
const FormFactors ela = FormFactors::ProtonElastic( Q2_ela );
const double ela1 = ( 1.-x )*( 1.-Q2_min/Q2_ela );
const double ela2 = ela.FE;
//const double ela3 = 1.-( Q2_ela-kt2 )/Q2_ela;
const double f_ela = Constants::alphaEM/M_PI/kt2*( ela1*ela2 + 0.5*x*x*ela.FM );
return f_ela * ( 1.-x ) * kt2 / Q2_ela;
}
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 ) + x*x*mp2 ),
Q2 = Q2min + kt2 / ( 1.-x );
float xbj = Q2 / ( Q2 + mx2 - mp2 );
StructureFunctions sf;
switch ( cuts_.structure_functions ) {
case StructureFunctionsType::SzczurekUleshchenko: {
SF::SzczurekUleshchenko su;
sf = su( Q2, xbj ); } break;
case StructureFunctionsType::SuriYennie: {
SF::SuriYennie sy;
sf = sy( Q2, xbj ); } break;
- case StructureFunctionsType::Fiore: {
+ case StructureFunctionsType::FioreBrasse: {
SF::FioreBrasse fb;
sf = fb( Q2, xbj ); } break;
case StructureFunctionsType::ALLM91: {
SF::ALLM allm91( SF::ALLM::Parameterisation::allm91() );
sf = allm91( Q2, xbj ); } break;
case StructureFunctionsType::ALLM97: {
SF::ALLM allm97( SF::ALLM::Parameterisation::allm97() );
sf = allm97( Q2, xbj ); } break;
+ case StructureFunctionsType::ALLM_GD07p: {
+ SF::ALLM gd07p( SF::ALLM::Parameterisation::gd07p() );
+ sf = gd07p( Q2, xbj ); } break;
+ case StructureFunctionsType::ALLM_GD11p: {
+ SF::ALLM gd11p( SF::ALLM::Parameterisation::gd11p() );
+ sf = gd11p( Q2, xbj ); } break;
default: break; //FIXME
}
// Longitudinal/transverse virtual photon cross section R
// from Sibirtsev and Blunden (Phys Rev C 88,065202 (2013))
const double R = 0.014 * Q2 * ( exp( -0.07*Q2 ) + 41.*exp( -0.8*Q2 ) );
const double F1 = sf.F2 * ( 1.+4*xbj*xbj*mp2/Q2 )/( 1.+R )/( 2.*xbj );
const double ine1 = ( 1.-x )*( 1.-Q2min / Q2 );
const double f_D = sf.F2/( Q2 + mx2 - mp2 ) * ine1, f_C= 2.*F1/Q2;
const double f_ine = Constants::alphaEM/M_PI/kt2*( f_D + 0.5*x*x*f_C );
return f_ine * ( 1.-x ) * kt2 / Q2;
}
}
}
diff --git a/CepGen/StructureFunctions/StructureFunctions.cpp b/CepGen/StructureFunctions/StructureFunctions.cpp
index 3bc4182..b3d2c0e 100644
--- a/CepGen/StructureFunctions/StructureFunctions.cpp
+++ b/CepGen/StructureFunctions/StructureFunctions.cpp
@@ -1,30 +1,28 @@
#include "StructureFunctions.h"
namespace CepGen
{
std::ostream&
operator<<( std::ostream& os, const StructureFunctionsType& sf )
{
switch ( sf ) {
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";
+ case FioreBrasse: return os << "Fiore-Brasse";
case ALLM91: return os << "ALLM;91";
case ALLM97: return os << "ALLM;97";
case ALLM_HHT: return os << "ALLM;HHT";
case ALLM_HHT_FT: return os << "ALLM;HHT-FT";
default: return os;
}
}
std::ostream&
operator<<( std::ostream& os, const StructureFunctions& sf )
{
return os << "F2 = " << sf.F2 << ", FL = " << sf.FL;
}
}
diff --git a/CepGen/StructureFunctions/StructureFunctions.h b/CepGen/StructureFunctions/StructureFunctions.h
index c8de7f8..32d082c 100644
--- a/CepGen/StructureFunctions/StructureFunctions.h
+++ b/CepGen/StructureFunctions/StructureFunctions.h
@@ -1,40 +1,38 @@
#ifndef CepGen_StructureFunctions_StructureFunctions_h
#define CepGen_StructureFunctions_StructureFunctions_h
#include "CepGen/Event/Particle.h"
#include <ostream>
#include <vector>
#include <cmath>
namespace CepGen
{
/// Proton structure function to be used in the outgoing state description
/// \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,
- ALLM91, ALLM97, ALLM_HHT, ALLM_HHT_FT
+ FioreBrasse = 101,
+ ALLM91, ALLM97, ALLM_HHT, ALLM_HHT_FT, ALLM_GD07p, ALLM_GD11p
};
/// Human-readable format of a structure function type
std::ostream& operator<<( std::ostream& os, const StructureFunctionsType& sf );
class StructureFunctions
{
public:
StructureFunctions( double f2=0.0 ) : F2( f2 ), FL( 0.0 ), FM( 0.0 ) {}
double F2;
double F1;
double FL, FM;
};
/// Human-readable format of a structure function object
std::ostream& operator<<( std::ostream& os, const StructureFunctions& sf );
}
#endif
diff --git a/CepGen/StructureFunctions/SuriYennie.cpp b/CepGen/StructureFunctions/SuriYennie.cpp
index 8b4318a..f2c18e1 100644
--- a/CepGen/StructureFunctions/SuriYennie.cpp
+++ b/CepGen/StructureFunctions/SuriYennie.cpp
@@ -1,51 +1,51 @@
#include "SuriYennie.h"
namespace CepGen
{
namespace SF
{
SuriYennie::Parameterisation
SuriYennie::Parameterisation::standard()
{
Parameterisation p;
p.C1 = 0.86926;
p.C2 = 2.23422;
p.D1 = 0.12549;
p.rho2 = 0.585;
p.Cp = 0.96;
p.Bp = 0.63;
return p;
}
SuriYennie::Parameterisation
SuriYennie::Parameterisation::alternative()
{
Parameterisation p;
p.C1 = 0.6303;
p.C2 = 2.3049;
p.D1 = 0.04681;
p.rho2 = 1.05;
p.Cp = 1.23;
p.Bp = 0.61;
return p;
}
StructureFunctions
SuriYennie::operator()( double q2, double xbj ) const
{
const double mp = Particle::massFromPDGId( Particle::Proton ), mp2 = mp*mp;
const double mx2 = q2 * ( 1.-xbj )/xbj + mp2, // [GeV^2]
nu = 0.5 * ( q2 + mx2 - mp2 ) / mp; // [GeV]
const double dm2 = mx2-mp2, Xpr = q2/( q2+mx2 ), En = dm2+q2, Tau = 0.25 * q2/mp2, MQ = params_.rho2+q2;
StructureFunctions sy;
sy.FM = ( params_.C1*dm2*pow( params_.rho2/MQ, 2 ) + ( params_.C2*mp2*pow( 1.-Xpr, 4 ) ) / ( 1.+Xpr*( Xpr*params_.Cp-2.*params_.Bp ) ) )/q2;
const double FE = ( Tau*sy.FM + params_.D1*dm2*q2*params_.rho2/mp2*pow( dm2/MQ/En, 2 ) ) / ( 1.+0.25*En*En/mp2/q2 );
- const double w2 = 2.*mp*FE, w1 = 0.5 * sy.FM*q2/mp;
+ const double w2 = 2.*mp*FE/*, w1 = 0.5 * sy.FM*q2/mp*/;
sy.F2 = nu/mp*w2;
return sy;
}
}
}