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( pname ) ) + return kParameters.get( pname ).pdgid; + return kParameters.get( pname, def ); + } + + double + cepgen_param_real_( char* pname, double& def ) + { + return kParameters.get( 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 + FortranKTProcess::kProcParameters; + FortranKTProcess::FortranKTProcess( const ParametersList& params, const char* name, const char* descr, std::function func ) : GenericKTProcess( params, name, descr, { { PDG::photon, PDG::photon } }, { PDG::muon, PDG::muon } ), - pair_( params.get( "pair" ).pdgid ), - method_( params.get( "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 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 func ); ProcessPtr clone( const ParametersList& /*params*/ ) const override { return ProcessPtr( new FortranKTProcess( *this ) ); } + static std::unordered_map 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 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( 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 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