diff --git a/CepGen/Core/FortranInterface.cpp b/CepGen/Core/FortranInterface.cpp
index 47346ee..f926d7c 100644
--- a/CepGen/Core/FortranInterface.cpp
+++ b/CepGen/Core/FortranInterface.cpp
@@ -1,76 +1,110 @@
#include "CepGen/StructureFunctions/StructureFunctions.h"
+#include "CepGen/Processes/FortranKTProcess.h"
#include "CepGen/Physics/KTFlux.h"
#include "CepGen/Physics/HeavyIon.h"
#include "CepGen/Physics/PDG.h"
#include "CepGen/Core/ParametersList.h"
#include "CepGen/Core/Exception.h"
#ifdef __cplusplus
extern "C" {
#endif
/// Expose structure functions calculators to Fortran
void
cepgen_structure_functions_( int& sfmode, double& xbj, double& q2, double& f2, double& fl )
{
using namespace cepgen;
static auto sf = strfun::StructureFunctionsHandler::get().build( sfmode );
const auto& val = ( *sf )( xbj, q2 );
f2 = val.F2;
fl = val.FL;
}
/// Compute a \f$k_{\rm T}\f$-dependent flux for single nucleons
/// \param[in] fmode Flux mode (see cepgen::KTFlux)
/// \param[in] x Fractional momentum loss
/// \param[in] kt2 The \f$k_{\rm T}\f$ transverse momentum norm
/// \param[in] sfmode Structure functions set for dissociative emission
/// \param[in] mx Diffractive state mass for dissociative emission
double
cepgen_kt_flux_( int& fmode, double& x, double& kt2, int& sfmode, double& mx )
{
using namespace cepgen;
static auto sf = strfun::StructureFunctionsHandler::get().build( sfmode );
return ktFlux( (KTFlux)fmode, x, kt2, *sf, mx );
}
/// Compute a \f$k_{\rm T}\f$-dependent flux for heavy ions
/// \param[in] fmode Flux mode (see cepgen::KTFlux)
/// \param[in] x Fractional momentum loss
/// \param[in] kt2 The \f$k_{\rm T}\f$ transverse momentum norm
/// \param[in] a Mass number for the heavy ion
/// \param[in] z Atomic number for the heavy ion
double
cepgen_kt_flux_hi_( int& fmode, double& x, double& kt2, int& a, int& z )
{
using namespace cepgen;
return ktFlux( (KTFlux)fmode, x, kt2, HeavyIon{ (unsigned short)a, (Element)z } );
}
/// Mass of a particle, in GeV/c^2
double
cepgen_particle_mass_( int& pdg_id )
{
try {
return cepgen::PDG::get().mass( (cepgen::pdgid_t)pdg_id );
} catch ( const cepgen::Exception& e ) {
e.dump();
exit( 0 );
}
}
/// Charge of a particle, in e
double
cepgen_particle_charge_( int& pdg_id )
{
try {
return cepgen::PDG::get().charge( (cepgen::pdgid_t)pdg_id );
} catch ( const cepgen::Exception& e ) {
e.dump();
exit( 0 );
}
}
+
+ static cepgen::ParametersList kParameters;
+
+ void
+ cepgen_set_process_( char* prname, int prlen )
+ {
+ const std::string proc( prname, prlen );
+ if ( cepgen::proc::FortranKTProcess::kProcParameters.count( proc ) == 0 )
+ throw CG_FATAL( "cepgen_set_process" )
+ << "Failed to retrieve a process named \"" << proc << "\" from the database!";
+ kParameters = cepgen::proc::FortranKTProcess::kProcParameters.at( proc );
+ }
+
+ void
+ cepgen_list_params_()
+ {
+ CG_LOG( "cepgen_list_params" ) << "\t" << kParameters;
+ }
+
+ int
+ cepgen_param_int_( char* pname, int& def )
+ {
+ if ( kParameters.has<cepgen::ParticleProperties>( pname ) )
+ return kParameters.get<cepgen::ParticleProperties>( pname ).pdgid;
+ return kParameters.get<int>( pname, def );
+ }
+
+ double
+ cepgen_param_real_( char* pname, double& def )
+ {
+ return kParameters.get<double>( pname, def );
+ }
+
#ifdef __cplusplus
}
#endif
diff --git a/CepGen/Processes/CMakeLists.txt b/CepGen/Processes/CMakeLists.txt
index 7489f0c..aadc93d 100644
--- a/CepGen/Processes/CMakeLists.txt
+++ b/CepGen/Processes/CMakeLists.txt
@@ -1,14 +1,15 @@
file(GLOB sources *.cpp Fortran/*.f)
include_directories(${PROJECT_SOURCE_DIR})
include_directories(Fortran)
set(F77_PROC_PATH ${PROJECT_SOURCE_DIR}/External/Processes)
file(GLOB oth_proc_sources ${F77_PROC_PATH}/*.f ${F77_PROC_PATH}/CepGenWrapper.cpp)
if(oth_proc_sources)
+ message(STATUS "Fortran processes found: ${oth_proc_sources}")
list(APPEND sources ${oth_proc_sources})
endif()
add_library(CepGenProcesses SHARED ${sources})
install(TARGETS CepGenProcesses DESTINATION lib)
diff --git a/CepGen/Processes/Fortran/KTBlocks.inc b/CepGen/Processes/Fortran/KTBlocks.inc
index da5c605..6b9a1b6 100644
--- a/CepGen/Processes/Fortran/KTBlocks.inc
+++ b/CepGen/Processes/Fortran/KTBlocks.inc
@@ -1,61 +1,68 @@
!> \file KTBblocks.inc
c ================================================================
c F77 process-to-CepGen helper
c > this header file defines a set of common blocks useful for
c > the interfacing of any kt-factorised matrix element computation
c > to a mother CepGen instance
c ================================================================
c =================================================================
!> collection of fundamental physics constants
common/constants/am_p,units,pi,alpha_em
double precision am_p,units,pi,alpha_em
c =================================================================
c information on the full process
c inp1 = proton energy in lab frame
c inp2 = nucleus energy **per nucleon** in LAB frame
c Collision is along z-axis
c =================================================================
- common/genparams/icontri,iflux1,iflux2,imethod,sfmod,pdg_l,
- & a_nuc1,z_nuc1,a_nuc2,z_nuc2,
- & inp1,inp2
- integer icontri,iflux1,iflux2,imethod,sfmod,pdg_l,
+ common/genparams/icontri,iflux1,iflux2,sfmod,
+ & a_nuc1,z_nuc1,a_nuc2,z_nuc2,inp1,inp2
+ integer icontri,iflux1,iflux2,sfmod,
& a_nuc1,z_nuc1,a_nuc2,z_nuc2
double precision inp1,inp2
c =================================================================
c kt-factorisation kinematics
c =================================================================
common/ktkin/q1t,q2t,phiq1t,phiq2t,y1,y2,ptdiff,phiptdiff,
& am_x,am_y
double precision q1t,q2t,phiq1t,phiq2t,y1,y2,ptdiff,phiptdiff,
& am_x,am_y
c =================================================================
c phase space cuts
c =================================================================
common/kincuts/ipt,iene,ieta,iinvm,iptsum,idely,
& pt_min,pt_max,ene_min,ene_max,eta_min,eta_max,
& invm_min,invm_max,ptsum_min,ptsum_max,
& dely_min,dely_max
logical ipt,iene,ieta,iinvm,iptsum,idely
double precision pt_min,pt_max,ene_min,ene_max,eta_min,eta_max,
& invm_min,invm_max,ptsum_min,ptsum_max,
& dely_min,dely_max
c =================================================================
c generated event kinematics
c =================================================================
common/evtkin/nout,ipdg,idum,pc,px,py
integer nout,idum,ipdg(10)
double precision pc(10,4),px(4),py(4)
c =================================================================
c helpers for the evaluation of unintegrated parton fluxes
c =================================================================
external CepGen_kT_flux,CepGen_kT_flux_HI
double precision CepGen_kT_flux,CepGen_kT_flux_HI
external CepGen_particle_charge,CepGen_particle_mass
double precision CepGen_particle_charge,CepGen_particle_mass
+c =================================================================
+c helpers for the retrieval of input parameters from the process
+c =================================================================
+ external CepGen_set_process
+ external CepGen_param_int,CepGen_param_real
+ integer CepGen_param_int
+ double precision CepGen_param_real
+
diff --git a/CepGen/Processes/Fortran/KTStructures.h b/CepGen/Processes/Fortran/KTStructures.h
index cad251b..a9d7b3f 100644
--- a/CepGen/Processes/Fortran/KTStructures.h
+++ b/CepGen/Processes/Fortran/KTStructures.h
@@ -1,77 +1,75 @@
#ifndef CepGen_Processes_Fortran_KTStructures_h
#define CepGen_Processes_Fortran_KTStructures_h
namespace cepgen
{
/// Collection of common blocks for Fortran \f$k_{\rm T}\f$-processes
namespace ktblock
{
/// General physics constants
struct Constants {
double m_p; ///< Proton mass
double units; ///< Conversion factor GeV\f$^2\to\f$ barn
double pi; ///< \f$\pi\f$
double alpha_em; ///< Electromagnetic coupling constant
};
/// Generic run parameters
struct Parameters {
int icontri; ///< Kinematics mode
int iflux1; ///< Type of \f$k_{\rm T}\f$-factorised flux for first incoming parton
int iflux2; ///< Type of \f$k_{\rm T}\f$-factorised flux for second incoming parton
- int imethod; ///< Computation method for matrix element
int sfmod; ///< Structure functions modelling
- int pdg_l; ///< Central system PDG id
int a_nuc1; ///< First beam mass number
int z_nuc1; ///< First beam atomic number
int a_nuc2; ///< Second beam mass number
int z_nuc2; ///< Second beam atomic number
double inp1; ///< First beam momentum, in GeV/c
double inp2; ///< Second beam momentum, in GeV/c
};
/// Kinematics properties of the \f$k_{\rm T}\f$-factorised process
struct KTKinematics {
double q1t; ///< Transverse momentum of the first incoming parton
double q2t; ///< Transverse momentum of the second incoming parton
double phiq1t; ///< Azimutal angle of the first incoming parton
double phiq2t; ///< Azimutal angle of the second incoming parton
double y1; ///< First incoming parton rapidity
double y2; ///< Second incoming parton rapidity
double ptdiff; ///< Central system pT balance
double phiptdiff; ///< Central system azimutal angle difference
double m_x; ///< Invariant mass for the first diffractive state
double m_y; ///< Invariant mass for the second diffractive state
};
/// Phase space cuts for event kinematics
struct Cuts {
bool ipt; ///< Switch for cut on single particle transverse momentum
bool iene; ///< Switch for cut on single particle energy
bool ieta; ///< Switch for cut on single particle pseudo-rapidity
bool iinvm; ///< Switch for cut on central system invariant mass
bool iptsum; ///< Switch for cut on central system transverse momentum
bool idely; ///< Switch for cut on rapididty difference
double pt_min; ///< Minimal single particle transverse momentum
double pt_max; ///< Maximal single particle transverse momentum
double ene_min; ///< Minimal single particle energy
double ene_max; ///< Maximal single particle energy
double eta_min; ///< Minimal single particle pseudo-rapidity
double eta_max; ///< Maximal single particle pseudo-rapidity
double invm_min; ///< Minimal central system invariant mass
double invm_max; ///< Maximal central system invariant mass
double ptsum_min; ///< Minimal central system transverse momentum
double ptsum_max; ///< Maximal central system transverse momentum
double dely_min; ///< Minimal rapidity difference for central system
double dely_max; ///< Maximal rapidity difference for central system
};
/// Single event kinematics
struct Event {
int nout; ///< Number of particles in central system
int pdg[10]; ///< PDG ids of all particles in central system
int idum; ///< Padding
double pc[4][10]; ///< 4-momenta of all particles in central system
double px[4]; ///< 4-momentum of first outgoing proton state
double py[4]; ///< 4-momentum of second outgoing proton state
};
}
}
#endif
diff --git a/CepGen/Processes/Fortran/cepgen_print.f b/CepGen/Processes/Fortran/cepgen_print.f
index af7a984..b99fe70 100644
--- a/CepGen/Processes/Fortran/cepgen_print.f
+++ b/CepGen/Processes/Fortran/cepgen_print.f
@@ -1,42 +1,43 @@
!> \file cepgen_print.f
subroutine cepgen_print
!> Print useful run information in standard stream
implicit none
include 'KTBlocks.inc'
logical params_shown
data params_shown/.false./
save params_shown
if(params_shown) return
print *,'========================================================'
print *,'Parameter value(s)'
print *,'--------------------------------------------------------'
print 101,'Process mode:',icontri
- print 101,'Computation method:',imethod
print 101,'Structure functions:',sfmod
- print 101,'Central system PDG:',pdg_l
print 103,'Beams momenta:',inp1,inp2
print 102,'Fluxes modes:',iflux1,iflux2
print 104,'Beams (A,Z):',a_nuc1,z_nuc1,a_nuc2,z_nuc2
print *,'========================================================'
print *,'Cut enabled minimum maximum'
print *,'--------------------------------------------------------'
print 100,'pt(single)',ipt,pt_min,pt_max
print 100,'energy(single)',iene,ene_min,ene_max
print 100,'eta(single)',ieta,eta_min,eta_max
print 100,'m(sum)',iinvm,invm_min,invm_max
print 100,'pt(sum)',iptsum,ptsum_min,ptsum_max
print 100,'delta(y)',idely,dely_min,dely_max
print *,'========================================================'
+ print *,'Process-specific parameters'
+ call cepgen_list_params
+ print *,'========================================================'
params_shown=.true.
100 format(A26,' ',L2,f12.4,f12.4)
101 format(A33,I12)
102 format(A33,I12,I12)
103 format(A33,f12.2,f12.2)
104 format(A33,' (',I3,',',I3,'), (',I3,',',I3,')')
end
diff --git a/CepGen/Processes/FortranKTProcess.cpp b/CepGen/Processes/FortranKTProcess.cpp
index aae8742..46e62ed 100644
--- a/CepGen/Processes/FortranKTProcess.cpp
+++ b/CepGen/Processes/FortranKTProcess.cpp
@@ -1,181 +1,180 @@
#include "CepGen/Processes/FortranKTProcess.h"
#include "CepGen/Processes/Fortran/KTStructures.h"
#include "CepGen/Core/ParametersList.h"
#include "CepGen/Core/Exception.h"
#include "CepGen/StructureFunctions/StructureFunctions.h"
#include "CepGen/Event/Event.h"
#include "CepGen/Physics/KTFlux.h"
#include "CepGen/Physics/Constants.h"
#include "CepGen/Physics/PDG.h"
extern "C"
{
extern cepgen::ktblock::Constants constants_;
extern cepgen::ktblock::Parameters genparams_;
extern cepgen::ktblock::KTKinematics ktkin_;
extern cepgen::ktblock::Cuts kincuts_;
extern cepgen::ktblock::Event evtkin_;
}
namespace cepgen
{
namespace proc
{
+ std::unordered_map<std::string,ParametersList>
+ FortranKTProcess::kProcParameters;
+
FortranKTProcess::FortranKTProcess( const ParametersList& params, const char* name, const char* descr, std::function<double( void )> func ) :
GenericKTProcess( params, name, descr, { { PDG::photon, PDG::photon } }, { PDG::muon, PDG::muon } ),
- pair_( params.get<ParticleProperties>( "pair" ).pdgid ),
- method_( params.get<int>( "method", 1 ) ),
func_( func )
{
constants_.m_p = GenericProcess::mp_;
constants_.units = constants::GEVM2_TO_PB;
constants_.pi = M_PI;
constants_.alpha_em = constants::ALPHA_EM;
}
void
FortranKTProcess::preparePhaseSpace()
{
mom_ip1_ = event_->getOneByRole( Particle::IncomingBeam1 ).momentum();
mom_ip2_ = event_->getOneByRole( Particle::IncomingBeam2 ).momentum();
registerVariable( y1_, Mapping::linear, kin_.cuts.central.rapidity_single, { -6., 6. }, "First central particle rapidity" );
registerVariable( y2_, Mapping::linear, kin_.cuts.central.rapidity_single, { -6., 6. }, "Second central particle rapidity" );
registerVariable( pt_diff_, Mapping::linear, kin_.cuts.central.pt_diff, { 0., 50. }, "Transverse momentum difference between central particles" );
registerVariable( phi_pt_diff_, Mapping::linear, kin_.cuts.central.phi_pt_diff, { 0., 2.*M_PI }, "Central particles azimuthal angle difference" );
//===========================================================================================
// feed phase space cuts to the common block
//===========================================================================================
kin_.cuts.central.pt_single.save( kincuts_.ipt, kincuts_.pt_min, kincuts_.pt_max );
kin_.cuts.central.energy_single.save( kincuts_.iene, kincuts_.ene_min, kincuts_.ene_max );
kin_.cuts.central.eta_single.save( kincuts_.ieta, kincuts_.eta_min, kincuts_.eta_max );
kin_.cuts.central.mass_sum.save( kincuts_.iinvm, kincuts_.invm_min, kincuts_.invm_max );
kin_.cuts.central.pt_sum.save( kincuts_.iptsum, kincuts_.ptsum_min, kincuts_.ptsum_max );
kin_.cuts.central.rapidity_diff.save( kincuts_.idely, kincuts_.dely_min, kincuts_.dely_max );
//===========================================================================================
// feed run parameters to the common block
//===========================================================================================
genparams_.icontri = (int)kin_.mode;
- genparams_.imethod = method_;
genparams_.sfmod = (int)kin_.structure_functions->type;
- genparams_.pdg_l = pair_;
//-------------------------------------------------------------------------------------------
// incoming beams information
//-------------------------------------------------------------------------------------------
genparams_.inp1 = kin_.incoming_beams.first.pz;
genparams_.inp2 = kin_.incoming_beams.second.pz;
const HeavyIon in1 = (HeavyIon)kin_.incoming_beams.first.pdg;
if ( in1 ) {
genparams_.a_nuc1 = in1.A;
genparams_.z_nuc1 = (unsigned short)in1.Z;
if ( genparams_.z_nuc1 > 1 ) {
event_->getOneByRole( Particle::IncomingBeam1 ).setPdgId( (pdgid_t)in1 );
event_->getOneByRole( Particle::OutgoingBeam1 ).setPdgId( (pdgid_t)in1 );
}
}
else
genparams_.a_nuc1 = genparams_.z_nuc1 = 1;
const HeavyIon in2 = (HeavyIon)kin_.incoming_beams.second.pdg;
if ( in2 ) {
genparams_.a_nuc2 = in2.A;
genparams_.z_nuc2 = (unsigned short)in2.Z;
if ( genparams_.z_nuc2 > 1 ) {
event_->getOneByRole( Particle::IncomingBeam2 ).setPdgId( (pdgid_t)in2 );
event_->getOneByRole( Particle::OutgoingBeam2 ).setPdgId( (pdgid_t)in2 );
}
}
else
genparams_.a_nuc2 = genparams_.z_nuc2 = 1;
//-------------------------------------------------------------------------------------------
// intermediate partons information
//-------------------------------------------------------------------------------------------
genparams_.iflux1 = (int)kin_.incoming_beams.first.kt_flux;
genparams_.iflux2 = (int)kin_.incoming_beams.second.kt_flux;
switch ( (KTFlux)genparams_.iflux1 ) {
case KTFlux::P_Gluon_KMR:
event_->getOneByRole( Particle::Parton1 ).setPdgId( PDG::gluon ); break;
case KTFlux::P_Photon_Elastic: case KTFlux::P_Photon_Inelastic: case KTFlux::P_Photon_Inelastic_Budnev:
case KTFlux::HI_Photon_Elastic:
event_->getOneByRole( Particle::Parton2 ).setPdgId( PDG::photon ); break;
case KTFlux::invalid:
throw CG_FATAL( "FortranKTProcess" )
<< "Invalid flux for 2nd incoming parton: " << genparams_.iflux2 << "!";
}
switch ( (KTFlux)genparams_.iflux2 ) {
case KTFlux::P_Gluon_KMR:
event_->getOneByRole( Particle::Parton2 ).setPdgId( PDG::gluon ); break;
case KTFlux::P_Photon_Elastic: case KTFlux::P_Photon_Inelastic: case KTFlux::P_Photon_Inelastic_Budnev:
case KTFlux::HI_Photon_Elastic:
event_->getOneByRole( Particle::Parton2 ).setPdgId( PDG::photon ); break;
case KTFlux::invalid:
throw CG_FATAL( "FortranKTProcess" )
<< "Invalid flux for 2nd incoming parton: " << genparams_.iflux2 << "!";
}
}
double
FortranKTProcess::computeKTFactorisedMatrixElement()
{
ktkin_.q1t = qt1_;
ktkin_.q2t = qt2_;
ktkin_.phiq1t = phi_qt1_;
ktkin_.phiq2t = phi_qt2_;
ktkin_.y1 = y1_;
ktkin_.y2 = y2_;
ktkin_.ptdiff = pt_diff_;
ktkin_.phiptdiff = phi_pt_diff_;
ktkin_.m_x = MX_;
ktkin_.m_y = MY_;
//--- compute the event weight
return func_();
}
void
FortranKTProcess::fillCentralParticlesKinematics()
{
//===========================================================================================
// outgoing beam remnants
//===========================================================================================
PX_ = Particle::Momentum( evtkin_.px );
PY_ = Particle::Momentum( evtkin_.py );
// express these momenta per nucleon
PX_ *= 1./genparams_.a_nuc1;
PY_ *= 1./genparams_.a_nuc2;
//===========================================================================================
// intermediate partons
//===========================================================================================
const Particle::Momentum mom_par1 = mom_ip1_-PX_, mom_par2 = mom_ip2_-PY_;
event_->getOneByRole( Particle::Parton1 ).setMomentum( mom_par1 );
event_->getOneByRole( Particle::Parton2 ).setMomentum( mom_par2 );
event_->getOneByRole( Particle::Intermediate ).setMomentum( mom_par1+mom_par2 );
//===========================================================================================
// central system
//===========================================================================================
Particles& oc = (*event_)[Particle::CentralSystem];
for ( int i = 0; i < evtkin_.nout; ++i ) {
Particle& p = oc[i];
p.setPdgId( (pdgid_t)evtkin_.pdg[i] );
p.setStatus( Particle::Status::FinalState );
p.setMomentum( Particle::Momentum( evtkin_.pc[i] ) );
}
}
}
}
diff --git a/CepGen/Processes/FortranKTProcess.h b/CepGen/Processes/FortranKTProcess.h
index ea6aa42..68a8519 100644
--- a/CepGen/Processes/FortranKTProcess.h
+++ b/CepGen/Processes/FortranKTProcess.h
@@ -1,37 +1,37 @@
#ifndef CepGen_Processes_FortranKTProcess_h
#define CepGen_Processes_FortranKTProcess_h
#include "CepGen/Processes/GenericKTProcess.h"
#include <functional>
namespace cepgen
{
namespace proc
{
/// Compute the matrix element for a generic \f$k_{\rm T}\f$-factorised process defined in a Fortran weighting function
class FortranKTProcess : public GenericKTProcess
{
public:
FortranKTProcess( const ParametersList& params, const char* name, const char* descr, std::function<double(void)> func );
ProcessPtr clone( const ParametersList& /*params*/ ) const override { return ProcessPtr( new FortranKTProcess( *this ) ); }
+ static std::unordered_map<std::string,ParametersList> kProcParameters;
+
private:
void preparePhaseSpace() override;
double computeKTFactorisedMatrixElement() override;
void fillCentralParticlesKinematics() override;
- int pair_; ///< Outgoing particles type
- int method_; ///< Computation method for the process
std::function<double(void)> func_; ///< Function to be called for weight computation
double y1_; ///< First outgoing particle rapidity
double y2_; ///< Second outgoing particle rapidity
double pt_diff_; ///< Transverse momentum balance between outgoing particles
double phi_pt_diff_; ///< Azimutal angle difference between outgoing particles
Particle::Momentum mom_ip1_; ///< First incoming beam momentum
Particle::Momentum mom_ip2_; ///< Second incoming beam momentum
};
}
}
#endif
diff --git a/CepGen/Processes/ProcessesHandler.h b/CepGen/Processes/ProcessesHandler.h
index 39603a7..93f03bc 100644
--- a/CepGen/Processes/ProcessesHandler.h
+++ b/CepGen/Processes/ProcessesHandler.h
@@ -1,40 +1,42 @@
#ifndef CepGen_Processes_ProcessesHandler_h
#define CepGen_Processes_ProcessesHandler_h
#include "CepGen/Core/ModuleFactory.h"
#include "CepGen/Core/ParametersList.h"
#include "CepGen/Processes/GenericProcess.h"
#include "CepGen/Processes/FortranKTProcess.h"
/** \file */
/// Add a generic process definition to the list of handled processes
#define REGISTER_PROCESS( name, obj ) \
namespace cepgen { namespace proc { \
struct BUILDERNM( name ) { \
BUILDERNM( name )() { ProcessesHandler::get().registerModule<obj>( STRINGIFY( name ) ); } }; \
static BUILDERNM( name ) g ## name; \
} }
/// Declare a Fortran process function name
#define DECLARE_FORTRAN_FUNCTION( method ) \
extern "C" { extern double method ## _(); }
#define PROCESS_F77_NAME( name ) F77_ ## name
/// Add the Fortran process definition to the list of handled processes
#define REGISTER_FORTRAN_PROCESS( name, method, description ) \
struct PROCESS_F77_NAME( name ) : public cepgen::proc::FortranKTProcess { \
PROCESS_F77_NAME( name )( const cepgen::ParametersList& params = cepgen::ParametersList() ) : \
- cepgen::proc::FortranKTProcess( params, STRINGIFY( name ), description, method ## _ ) {} }; \
+ cepgen::proc::FortranKTProcess( params, STRINGIFY( name ), description, method ## _ ) { \
+ cepgen::proc::FortranKTProcess::kProcParameters[STRINGIFY( method )] = params; \
+ } }; \
REGISTER_PROCESS( name, PROCESS_F77_NAME( name ) )
namespace cepgen
{
namespace proc
{
/// A processes factory
typedef ModuleFactory<GenericProcess> ProcessesHandler;
}
}
#endif
diff --git a/External/Processes/nucl_to_ff.f b/External/Processes/nucl_to_ff.f
index f5b65ce..b37f636 100644
--- a/External/Processes/nucl_to_ff.f
+++ b/External/Processes/nucl_to_ff.f
@@ -1,495 +1,498 @@
function nucl_to_ff()
implicit none
double precision nucl_to_ff
c =================================================================
c CepGen common blocks for kinematics definition
c =================================================================
include 'KTBlocks.inc'
data iflux1,iflux2,sfmod,pdg_l/10,100,11,13/
data a_nuc1,z_nuc1,a_nuc2,z_nuc2/1,1,208,82/
c =================================================================
c local variables
c =================================================================
double precision s,s12
double precision alpha1,alpha2,amt1,amt2
double precision pt1,pt2,eta1,eta2,dely
double precision pt1x,pt1y,pt2x,pt2y
double precision ak10,ak1x,ak1y,ak1z,ak20,ak2x,ak2y,ak2z
double precision beta1,beta2,x1,x2
double precision z1p,z1m,z2p,z2m
double precision q1tx,q1ty,q1z,q10,q1t2
double precision q2tx,q2ty,q2z,q20,q2t2
double precision ptsum
double precision that1,that2,that,uhat1,uhat2,uhat
double precision term1,term2,term3,term4,term5,term6,term7
double precision term8,term9,term10,amat2
double precision ak1_x,ak1_y,ak2_x,ak2_y
double precision eps12,eps22
double precision aux2_1,aux2_2,f1,f2
double precision Phi10,Phi102,Phi11_x,Phi11_y,Phi112
double precision Phi20,Phi202,Phi21_x,Phi21_y,Phi212
double precision Phi11_dot_e,Phi11_cross_e
double precision Phi21_dot_e,Phi21_cross_e
double precision aintegral
- integer imat1,imat2
+ integer imethod,pdg_l,imat1,imat2
double precision px_plus,px_minus,py_plus,py_minus
double precision r1,r2
double precision ptdiffx,ptsumx,ptdiffy,ptsumy
double precision invm,invm2,s1_eff,s2_eff
double precision t1abs,t2abs
double precision am_l,q_l
double precision inp_A,inp_B,am_A,am_B
double precision p1_plus, p2_minus
double precision amat2_1,amat2_2
integer iterm11,iterm22,iterm12,itermtt
double precision coupling
c =================================================================
c quarks production
c =================================================================
#ifdef ALPHA_S
double precision t_max,amu2,alphas
#endif
logical first_init
data first_init/.true./
save first_init,am_l,q_l
+ call CepGen_set_process('nucl_to_ff')
call CepGen_print
+ imethod = CepGen_param_int('method', 1)
+ pdg_l = CepGen_param_int('pair', 13)
if(first_init) then
am_l = CepGen_particle_mass(pdg_l) ! central particles mass
q_l = CepGen_particle_charge(pdg_l) ! central particles charge
if(iflux1.ge.20.and.iflux1.lt.40) then
if(icontri.eq.3.or.icontri.eq.4) then
print *,'Invalid process mode for collinear gluon emission!'
stop
endif
#ifdef ALPHA_S
print *,'Initialisation of the alpha(S) evolution algorithm..'
call initAlphaS(0,1.d0,1.d0,0.5d0,
& CepGen_particle_mass(4), ! charm
& CepGen_particle_mass(5), ! bottom
& CepGen_particle_mass(6)) ! top
#endif
endif
first_init = .false.
endif
c =================================================================
c FIXME
c =================================================================
nucl_to_ff = 0.d0
amat2 = 0.d0
eps12 = 0.d0
eps22 = 0.d0
q10 = 0.d0
q1z = 0.d0
q20 = 0.d0
q2z = 0.d0
c =================================================================
c go to energy of Nucleus = A*inp in order that generator puts out
c proper momenta in LAB frame
c =================================================================
inp_A = inp1*a_nuc1
am_A = am_p*a_nuc1
inp_B = inp2*a_nuc2
am_B = am_p*a_nuc2
c =================================================================
c four-momenta for incoming beams in LAB !!!!
c =================================================================
r1 = dsqrt(1.d0+am_A**2/inp_A**2)
r2 = dsqrt(1.d0+am_B**2/inp_B**2)
ak10 = inp_A*r1
ak1x = 0.d0
ak1y = 0.d0
ak1z = inp_A
ak20 = inp_B*r2
ak2x = 0.d0
ak2y = 0.d0
ak2z = -inp_B
s = 4.*inp_A*inp_B*(1.d0 + r1*r2)/2.d0+(am_A**2+am_B**2)
s12 = dsqrt(s)
p1_plus = (ak10+ak1z)/dsqrt(2.d0)
p2_minus = (ak20-ak2z)/dsqrt(2.d0)
c terms in the matrix element
c
iterm11 = 1 ! LL
iterm22 = 1 ! TT
iterm12 = 1 ! LT
itermtt = 1 ! TT'
c =================================================================
c two terms in the Wolfgang's formula for
c off-shell gamma gamma --> l^+ l^-
c =================================================================
imat1 = 1
imat2 = 1
c =================================================================
c Outgoing proton final state's mass
c =================================================================
if((icontri.eq.1).or.(icontri.eq.2)) am_x = am_A
if((icontri.eq.1).or.(icontri.eq.3)) am_y = am_B
q1tx = q1t*cos(phiq1t)
q1ty = q1t*sin(phiq1t)
q2tx = q2t*cos(phiq2t)
q2ty = q2t*sin(phiq2t)
ptsumx = q1tx+q2tx
ptsumy = q1ty+q2ty
ptsum = sqrt(ptsumx**2+ptsumy**2)
c =================================================================
c a window in final state transverse momentum
c =================================================================
if(iptsum) then
if(ptsum.lt.ptsum_min.or.ptsum.gt.ptsum_max) return
endif
c =================================================================
c compute the individual central particles momentum
c =================================================================
ptdiffx = ptdiff*cos(phiptdiff)
ptdiffy = ptdiff*sin(phiptdiff)
pt1x = 0.5*(ptsumx+ptdiffx)
pt1y = 0.5*(ptsumy+ptdiffy)
pt2x = 0.5*(ptsumx-ptdiffx)
pt2y = 0.5*(ptsumy-ptdiffy)
pt1 = sqrt(pt1x**2+pt1y**2)
pt2 = sqrt(pt2x**2+pt2y**2)
if(ipt) then
if(pt1.lt.pt_min.or.pt2.lt.pt_min) return
if(pt_max.gt.0d0.and.(pt2.gt.pt_max.or.pt2.gt.pt_max)) return
endif
amt1 = dsqrt(pt1**2+am_l**2)
amt2 = dsqrt(pt2**2+am_l**2)
invm2 = amt1**2 + amt2**2 + 2.d0*amt1*amt2*dcosh(y1-y2)
& -ptsum**2
invm = dsqrt(invm2)
c =================================================================
c a window in final state invariant mass
c =================================================================
if(iinvm) then
if(invm.lt.invm_min) return
if(invm_max.gt.0d0.and.invm.gt.invm_max) return
endif
c =================================================================
c a window in rapidity distance
c =================================================================
dely = dabs(y1-y2)
if(idely) then
if(dely.lt.dely_min.or.dely.gt.dely_max) return
endif
c =================================================================
c auxiliary quantities
c =================================================================
alpha1 = amt1/(dsqrt(2.d0)*p1_plus)*dexp( y1)
alpha2 = amt2/(dsqrt(2.d0)*p1_plus)*dexp( y2)
beta1 = amt1/(dsqrt(2.d0)*p2_minus)*dexp(-y1)
beta2 = amt2/(dsqrt(2.d0)*p2_minus)*dexp(-y2)
q1t2 = q1tx**2 + q1ty**2
q2t2 = q2tx**2 + q2ty**2
x1 = alpha1 + alpha2
x2 = beta1 + beta2
z1p = alpha1/x1
z1m = alpha2/x1
z2p = beta1/x2
z2m = beta2/x2
c -----------------------------------------------------------------
if(x1.gt.1.0.or.x2.gt.1.0) return
c -----------------------------------------------------------------
s1_eff = x1*s - q1t**2
s2_eff = x2*s - q2t**2
c -----------------------------------------------------------------
c Additional conditions for energy-momentum conservation
c -----------------------------------------------------------------
if(((icontri.eq.2).or.(icontri.eq.4))
1 .and.(dsqrt(s1_eff).le.(am_y+invm))) return
if(((icontri.eq.3).or.(icontri.eq.4))
1 .and.(dsqrt(s2_eff).le.(am_x+invm))) return
c =================================================================
c >>> TO THE OUTPUT COMMON BLOCK
c =================================================================
c =================================================================
c four-momenta of the outgoing protons (or remnants)
c =================================================================
px_plus = (1.d0-x1) * p1_plus
px_minus = ((a_nuc1*am_x)**2 + q1tx**2 + q1ty**2)/2.d0/px_plus
px(1) = - q1tx
px(2) = - q1ty
px(3) = (px_plus - px_minus)/dsqrt(2.d0)
px(4) = (px_plus + px_minus)/dsqrt(2.d0)
py_minus = (1.d0-x2) * p2_minus
py_plus = ((a_nuc2*am_y)**2 + q2tx**2 + q2ty**2)/2.d0/py_minus
py(1) = - q2tx
py(2) = - q2ty
py(3) = (py_plus - py_minus)/dsqrt(2.d0)
py(4) = (py_plus + py_minus)/dsqrt(2.d0)
q1t = dsqrt(q1t2)
q2t = dsqrt(q2t2)
c =================================================================
c four-momenta of the outgoing central particles
c =================================================================
nout = 2
ipdg(1) = pdg_l
pc(1,1) = pt1x
pc(1,2) = pt1y
pc(1,3) = alpha1*ak1z + beta1*ak2z
pc(1,4) = alpha1*ak10 + beta1*ak20
ipdg(2) = -pdg_l
pc(2,1) = pt2x
pc(2,2) = pt2y
pc(2,3) = alpha2*ak1z + beta2*ak2z
pc(2,4) = alpha2*ak10 + beta2*ak20
eta1 = 0.5d0*dlog((dsqrt(amt1**2*(dcosh(y1))**2 - am_l**2) +
2 amt1*dsinh(y1))/(dsqrt(amt1**2*(dcosh(y1))**2 - am_l**2)
3 - amt1*dsinh(y1)))
eta2 = 0.5d0*dlog((dsqrt(amt2**2*(dcosh(y2))**2 - am_l**2) +
2 amt2*dsinh(y2))/(dsqrt(amt2**2*(dcosh(y2))**2 - am_l**2)
3 - amt2*dsinh(y2)))
if(ieta) then
if(eta1.lt.eta_min.or.eta1.gt.eta_max) return
if(eta2.lt.eta_min.or.eta2.gt.eta_max) return
endif
c =================================================================
c matrix element squared
c averaged over initial spin polarizations
c and summed over final spin polarizations
c (--> see Wolfgang's notes
c =================================================================
c =================================================================
c Mendelstam variables
c =================================================================
that1 = (q10-pc(1,4))**2
& -(q1tx-pc(1,1))**2-(q1ty-pc(2,1))**2-(q1z-pc(3,1))**2
uhat1 = (q10-pc(2,4))**2
& -(q1tx-pc(1,2))**2-(q1ty-pc(2,2))**2-(q1z-pc(3,2))**2
that2 = (q20-pc(2,4))**2
& -(q2tx-pc(1,2))**2-(q2ty-pc(2,2))**2-(q2z-pc(3,2))**2
uhat2 = (q20-pc(1,4))**2
& -(q2tx-pc(1,1))**2-(q2ty-pc(2,1))**2-(q2z-pc(3,1))**2
that = (that1+that2)/2.d0
uhat = (uhat1+uhat2)/2.d0
c =================================================================
c matrix elements
c =================================================================
c How matrix element is calculated
c imethod = 0: on-shell formula
c imethod = 1: off-shell formula
c =================================================================
if(imethod.eq.0) then
c =================================================================
c on-shell formula for M^2
c =================================================================
term1 = 6.d0*am_l**8
term2 = -3.d0*am_l**4*that**2
term3 = -14.d0*am_l**4*that*uhat
term4 = -3.d0*am_l**4*uhat**2
term5 = am_l**2*that**3
term6 = 7.d0*am_l**2*that**2*uhat
term7 = 7.d0*am_l**2*that*uhat**2
term8 = am_l**2*uhat**3
term9 = -that**3*uhat
term10 = -that*uhat**3
amat2 = -2.d0*( term1+term2+term3+term4+term5
2 +term6+term7+term8+term9+term10 )
3 / ( (am_l**2-that)**2 * (am_l**2-uhat)**2 )
elseif(imethod.eq.1)then
c =================================================================
c Wolfgang's formulae
c =================================================================
ak1_x = z1m*pt1x-z1p*pt2x
ak1_y = z1m*pt1y-z1p*pt2y
ak2_x = z2m*pt1x-z2p*pt2x
ak2_y = z2m*pt1y-z2p*pt2y
!FIXME FIXME FIXME am_p or am_A/B???
t1abs = (q1t2+x1*(am_x**2-am_p**2)+x1**2*am_p**2)/(1.d0-x1)
t2abs = (q2t2+x2*(am_y**2-am_p**2)+x2**2*am_p**2)/(1.d0-x2)
eps12 = am_l**2 + z1p*z1m*t1abs
eps22 = am_l**2 + z2p*z2m*t2abs
Phi10 = 1.d0/((ak1_x+z1p*q2tx)**2+(ak1_y+z1p*q2ty)**2+eps12)
2 - 1.d0/((ak1_x-z1m*q2tx)**2+(ak1_y-z1m*q2ty)**2+eps12)
Phi11_x = (ak1_x+z1p*q2tx)/
2 ((ak1_x+z1p*q2tx)**2+(ak1_y+z1p*q2ty)**2+eps12)
3 - (ak1_x-z1m*q2tx)/
4 ((ak1_x-z1m*q2tx)**2+(ak1_y-z1m*q2ty)**2+eps12)
Phi11_y = (ak1_y+z1p*q2ty)/
2 ((ak1_x+z1p*q2tx)**2+(ak1_y+z1p*q2ty)**2+eps12)
3 - (ak1_y-z1m*q2ty)/
4 ((ak1_x-z1m*q2tx)**2+(ak1_y-z1m*q2ty)**2+eps12)
Phi102 = Phi10*Phi10
Phi112 = Phi11_x**2+Phi11_y**2
Phi20 = 1.d0/((ak2_x+z2p*q1tx)**2+(ak2_y+z2p*q1ty)**2+eps22)
2 - 1.d0/((ak2_x-z2m*q1tx)**2+(ak2_y-z2m*q1ty)**2+eps22)
Phi21_x = (ak2_x+z2p*q1tx)/
2 ((ak2_x+z2p*q1tx)**2+(ak2_y+z2p*q1ty)**2+eps22)
3 - (ak2_x-z2m*q1tx)/
4 ((ak2_x-z2m*q1tx)**2+(ak2_y-z2m*q1ty)**2+eps22)
Phi21_y = (ak2_y+z2p*q1ty)/
2 ((ak2_x+z2p*q1tx)**2+(ak2_y+z2p*q1ty)**2+eps22)
3 - (ak2_y-z2m*q1ty)/
4 ((ak2_x-z2m*q1tx)**2+(ak2_y-z2m*q1ty)**2+eps22)
Phi202 = Phi20*Phi20
Phi212 = Phi21_x**2+Phi21_y**2
Phi11_dot_e = (Phi11_x*q1tx + Phi11_y*q1ty)/dsqrt(q1t2)
Phi11_cross_e = (Phi11_x*q1ty - Phi11_y*q1tx)/dsqrt(q1t2)
Phi21_dot_e = (Phi21_x*q2tx +Phi21_y*q2ty)/dsqrt(q2t2)
Phi21_cross_e = (Phi21_x*q2ty -Phi21_y*q2tx)/dsqrt(q2t2)
aux2_1 = iterm11*(am_l**2+4.d0*z1p**2*z1m**2*t1abs)*Phi102
1 +iterm22*( (z1p**2 + z1m**2)*(Phi11_dot_e**2 +
2 Phi11_cross_e**2) )
3 + itermtt*( Phi11_cross_e**2 - Phi11_dot_e**2)
4 - iterm12*4.d0*z1p*z1m*(z1p-z1m)*Phi10
5 *(q1tx*Phi11_x+q1ty*Phi11_y)
aux2_2 = iterm11*(am_l**2+4.d0*z2p**2*z2m**2*t2abs)*Phi202
1 +iterm22*( (z2p**2 + z2m**2)*(Phi21_dot_e**2 +
2 Phi21_cross_e**2) )
3 + itermtt*( Phi21_cross_e**2 - Phi21_dot_e**2)
4 - iterm12*4.d0*z2p*z2m*(z2p-z2m)*Phi20
5 *(q2tx*Phi21_x+q2ty*Phi21_y)
c =================================================================
c convention of matrix element as in our kt-factorization
c for heavy flavours
c =================================================================
amat2_1 = (x1*x2*s)**2 * aux2_1 * 2.*z1p*z1m*q1t2 / (q1t2*q2t2)
amat2_2 = (x1*x2*s)**2 * aux2_2 * 2.*z2p*z2m*q2t2 / (q1t2*q2t2)
c =================================================================
c symmetrization
c =================================================================
amat2 = (imat1*amat2_1 + imat2*amat2_2)/2.d0
endif
coupling = 1.d0
c =================================================================
c first parton coupling
c =================================================================
if(iflux1.ge.20.and.iflux1.lt.40) then ! at least one gluon exchanged
#ifdef ALPHA_S
t_max = max(amt1**2,amt2**2)
amu2 = max(eps12,t_max)
am_x = dsqrt(amu2)
coupling = coupling * 4.d0*pi*alphaS(am_x)/2.d0 ! colour flow
#else
print *,'alphaS not linked to this instance!'
stop
#endif
else ! photon exchange
coupling = coupling * 4.d0*pi*alpha_em*q_l**2
endif
c =================================================================
c second parton coupling
c =================================================================
coupling = coupling * 4.d0*pi*alpha_em*q_l**2 ! photon exchange
coupling = coupling * 3.d0
c ============================================
c unintegrated parton distributions
c ============================================
if(a_nuc1.le.1) then
f1 = CepGen_kT_flux(iflux1,x1,q1t2,sfmod,am_x)
else
f1 = CepGen_kT_flux_HI(iflux1,x1,q1t2,a_nuc1,z_nuc1)
endif
if(a_nuc2.le.1) then
f2 = CepGen_kT_flux(iflux2,x2,q2t2,sfmod,am_y)
else
f2 = CepGen_kT_flux_HI(iflux2,x2,q2t2,a_nuc2,z_nuc2)
endif
c =================================================================
c factor 2.*pi below from integration over phi_sum
c factor 1/4 below from jacobian of transformations
c factors 1/pi and 1/pi due to integration
c over d^2 kappa_1 d^2 kappa_2 instead d kappa_1^2 d kappa_2^2
c =================================================================
aintegral = (2.d0*pi)*1.d0/(16.d0*pi**2*(x1*x2*s)**2) * amat2
& * coupling
& * f1/pi * f2/pi * (1.d0/4.d0) * units
& * 0.5d0 * 4.0d0 / (4.d0*pi)
c *****************************************************************
c =================================================================
nucl_to_ff = aintegral*q1t*q2t*ptdiff
c =================================================================
c *****************************************************************
c print *,nucl_to_ff,aintegral,coupling
return
end