Page MenuHomeHEPForge

No OneTemporary

diff --git a/Cards/pptoll.cfg b/Cards/pptoll.cfg
index d64d6d5..0572aff 100644
--- a/Cards/pptoll.cfg
+++ b/Cards/pptoll.cfg
@@ -1,35 +1,35 @@
@include "Cards/vegas.cfi"
process = {
name = "pptoll";
mode = "elastic/elastic";
in_kinematics = {
beam1_pz = 6500.;
beam2_pz = 6500.;
structure_functions = "Suri-Yennie";
};
out_kinematics = {
pair = 13;
min_pt = 25.0;
min_energy = 0.0;
- min_eta = -2.5;
- max_eta = 2.5;
+ min_rapidity = -2.5;
+ max_rapidity = 2.5;
max_mx = 1000.;
//--- extra cuts on the p1t(l) and p2t(l) plane
//max_ptdiff = 2.5;
//--- distance in rapidity between l^+ and l^-
//min_dely = 4.0;
//max_dely = 5.0;
};
rescattering_corr = "min(1,exp(-x/10))";
};
//--- either use the default generation (100k events)
@include "Cards/generator.cfi"
//--- or let the user specify the run conditions
/*generator = {
num_events = 100000;
print_every = 10000;
};*/
diff --git a/CepGen/Cards/ConfigHandler.cpp b/CepGen/Cards/ConfigHandler.cpp
index 916b44e..ee82ad7 100644
--- a/CepGen/Cards/ConfigHandler.cpp
+++ b/CepGen/Cards/ConfigHandler.cpp
@@ -1,233 +1,235 @@
#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 std::string proc_name = proc["name"];
if ( proc_name == "lpair" ) params_.setProcess( new Process::GamGamLL );
else if ( proc_name == "pptoll" ) params_.setProcess( new Process::PPtoLL );
else FatalError( Form( "Unrecognised process: %s", proc_name.c_str() ) );
//--- process mode
int int_mode; std::string str_mode;
if ( proc.lookupValue( "mode", int_mode ) ) {
params_.kinematics.mode = (Kinematics::ProcessMode)int_mode;
}
else if ( proc.lookupValue( "mode", str_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", str_mode.c_str() ) );
}
//--- 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" ) ) {
std::string sf = kin["structure_functions" ];
if ( sf == "electron" ) params_.kinematics.structure_functions = Electron;
else if ( sf == "elastic proton" ) params_.kinematics.structure_functions = ElasticProton;
else if ( sf == "Suri-Yennie" ) params_.kinematics.structure_functions = SuriYennie;
else if ( sf == "Suri-Yennie;lowQ2" ) params_.kinematics.structure_functions = SuriYennieLowQ2;
else if ( sf == "Szczurek-Uleshchenko" ) params_.kinematics.structure_functions = SzczurekUleshchenko;
else if ( sf == "Fiore;valence" ) params_.kinematics.structure_functions = FioreVal;
else if ( sf == "Fiore;sea" ) params_.kinematics.structure_functions = FioreSea;
else if ( sf == "Fiore" ) params_.kinematics.structure_functions = Fiore;
else FatalError( Form( "Invalid structure functions mode: %s", sf.c_str() ) );
}
} 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"];
} 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 );
kin.add( "pair", libconfig::Setting::TypeInt ) = (int)params->kinematics.central_system[0];
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();
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();
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();
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();
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();
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;
}
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/Cuts.h b/CepGen/Physics/Cuts.h
index b423912..e62073f 100644
--- a/CepGen/Physics/Cuts.h
+++ b/CepGen/Physics/Cuts.h
@@ -1,67 +1,69 @@
#ifndef CepGen_Physics_Cuts_h
#define CepGen_Physics_Cuts_h
namespace CepGen
{
/// Constraints to be applied on the events kinematics
namespace Cuts
{
/// Cuts on the central particles (e.g. dilepton system for LPAIR)
enum Central
{
- pt_single, ///< single particle transverse momentum
- eta_single, ///< single particle pseudo-rapidity
- energy_single, ///< single particle energy
- mass_sum, ///< central system invariant mass
- pt_sum, ///< central system transverse momentum
- eta_sum, ///< central system pseudo-rapidity
- energy_sum, ///< central system energy
- pt_diff, ///< transverse momentum balance between the central particles
- dely ///< rapidity balance between the central particles
+ pt_single, ///< single particle transverse momentum
+ eta_single, ///< single particle pseudo-rapidity
+ rapidity_single, ///< single particle rapidity
+ energy_single, ///< single particle energy
+ mass_sum, ///< central system invariant mass
+ pt_sum, ///< central system transverse momentum
+ eta_sum, ///< central system pseudo-rapidity
+ energy_sum, ///< central system energy
+ pt_diff, ///< transverse momentum balance between the central particles
+ dely ///< rapidity balance between the central particles
};
inline std::ostream& operator<<( std::ostream& os, const Central& is ) {
switch ( is ) {
case pt_single: return os << "Single central particle pT";
case eta_single: return os << "Single central particle eta";
+ case rapidity_single: return os << "Single central particle rapidity";
case energy_single: return os << "Single central particle energy";
case mass_sum: return os << "Central system mass";
case pt_sum: return os << "Central system pT";
case eta_sum: return os << "Central system eta";
case energy_sum: return os << "Central system energy";
case pt_diff: return os << "Central system d(pT)";
case dely: return os << "Central system d(Y)";
}
return os;
}
/// Cuts on the beam particles remnants
enum Remnants
{
mass ///< outgoing beam particles remnants mass
};
inline std::ostream& operator<<( std::ostream& os, const Remnants& is ) {
switch ( is ) {
case mass: return os << "Remnants mass";
}
return os;
}
/// Cuts on the initial state dynamics
enum InitialState
{
q2, ///< parton virtuality
qt, ///< parton transverse virtuality
w ///< two-parton squared momentum
};
inline std::ostream& operator<<( std::ostream& os, const InitialState& is ) {
switch ( is ) {
case q2: return os << "Virtuality range";
case qt: return os << "Transverse virtuality range";
case w: return os << "W";
}
return os;
}
}
}
#endif
diff --git a/CepGen/Processes/GamGamLL.cpp b/CepGen/Processes/GamGamLL.cpp
index 4876847..0810831 100644
--- a/CepGen/Processes/GamGamLL.cpp
+++ b/CepGen/Processes/GamGamLL.cpp
@@ -1,877 +1,876 @@
#include "GamGamLL.h"
using namespace CepGen::Process;
GamGamLL::GamGamLL( int nopt ) :
GenericProcess( "lpair", "pp -> p(*) (gamma gamma -> l+ l-) p(*)" ),
n_opt_( nopt ),
MX2_( 0. ), MY2_( 0. ), Ml2_( 0. ),
ep1_( 0. ), ep2_( 0. ), p_cm_( 0. ),
w12_( 0. ), w31_( 0. ), dw31_( 0. ), w52_( 0. ), dw52_( 0. ),
ec4_( 0. ), pc4_( 0. ), mc4_( 0. ), w4_( 0. ),
p12_( 0. ), p1k2_( 0. ), p2k1_( 0. ),
p13_( 0. ), p14_( 0. ), p25_( 0. ),
q1dq_( 0. ), q1dq2_( 0. ),
s1_( 0. ), s2_( 0. ),
epsi_( 0. ),
g5_( 0. ), g6_( 0. ), a5_( 0. ), a6_( 0. ), bb_( 0. ),
gram_( 0. ),
dd1_( 0. ), dd2_( 0. ), dd3_( 0. ), dd4_( 0. ), dd5_( 0. ),
delta_( 0. ),
g4_( 0. ), sa1_( 0. ), sa2_( 0. ),
sl1_( 0. ),
cos_theta4_( 0. ), sin_theta4_( 0. ),
al4_( 0. ), be4_( 0. ), de3_( 0. ), de5_( 0. ),
pt4_( 0. ),
- jacobian_( 0. ),
- cot_theta1_( -99999. ), cot_theta2_( 99999. )
+ jacobian_( 0. )
{}
void
GamGamLL::addEventContent()
{
GenericProcess::setEventContent( {
{ Particle::IncomingBeam1, Particle::Proton },
{ Particle::IncomingBeam2, Particle::Proton },
{ Particle::Parton1, Particle::Photon },
{ Particle::Parton2, Particle::Photon }
}, {
{ Particle::OutgoingBeam1, { Particle::Proton } },
{ Particle::OutgoingBeam2, { Particle::Proton } },
{ Particle::CentralSystem, { Particle::Muon, Particle::Muon } }
} );
}
unsigned int
GamGamLL::numDimensions( const Kinematics::ProcessMode& process_mode ) const
{
switch ( process_mode ) {
case Kinematics::ElectronProton: { InError( "Not supported yet!" ); }
case Kinematics::ElasticElastic:
default: return 7;
case Kinematics::ElasticInelastic:
case Kinematics::InelasticElastic: return 8;
case Kinematics::InelasticInelastic: return 9;
}
}
bool
GamGamLL::pickin()
{
DebuggingInsideLoop( Form( "Optimised mode? %i", n_opt_ ) );
jacobian_ = 0.;
w4_ = mc4_*mc4_;
// sig1 = sigma and sig2 = sigma' in [1]
const double sig = mc4_+MY_;
double sig1 = sig*sig,
sig2 = sig1;
DebuggingInsideLoop( Form( "mc4 = %f\n\t"
"sig1 = %f\n\t"
"sig2 = %f", mc4_, sig1, sig2 ) );
// Mass difference between the first outgoing particle and the first incoming
// particle
w31_ = MX2_-w1_;
// Mass difference between the second outgoing particle and the second
// incoming particle
w52_ = MY2_-w2_;
// Mass difference between the two incoming particles
w12_ = w1_-w2_;
// Mass difference between the central two-photons system and the second
// outgoing particle
const double d6 = w4_-MY2_;
DebuggingInsideLoop( Form( "w1 = %f\n\t"
"w2 = %f\n\t"
"w3 = %f\n\t"
"w4 = %f\n\t"
"w5 = %f",
w1_, w2_, MX2_, w4_, MY2_ ) );
DebuggingInsideLoop( Form( "w31 = %f\n\tw52 = %f\n\tw12 = %f", w31_, w52_, w12_ ) );
const double ss = s_+w12_;
const double rl1 = ss*ss-4.*w1_*s_; // lambda(s, m1**2, m2**2)
if ( rl1 <= 0. ) { InWarning( Form( "rl1 = %f <= 0", rl1 ) ); return false; }
sl1_ = sqrt( rl1 );
s2_ = 0.;
double ds2 = 0.;
if ( n_opt_ == 0 ) {
const double smax = s_+MX2_-2.*MX_*sqs_;
Map( x(2), sig1, smax, s2_, ds2, "s2" );
sig1 = s2_; //FIXME!!!!!!!!!!!!!!!!!!!!
}
DebuggingInsideLoop( Form( "s2 = %f", s2_ ) );
//std::cout << "s=" << _s << ", w3=" << MX2_ << ", sig1=" << sig1 << std::endl;
const double sp = s_+MX2_-sig1,
d3 = sig1-w2_;
const double rl2 = sp*sp-4.*s_*MX2_; // lambda(s, m3**2, sigma)
if ( rl2 <= 0. ) { InWarning( Form( "rl2 = %f <= 0", rl2 ) ); return false; }
const double sl2 = sqrt( rl2 );
//std::cout << "ss=" << ss << ", sp=" << sp << ", sl1=" << sl1_ << ", sl2=" << sl2 << std::endl;
double t1_max = w1_+MX2_-( ss*sp+sl1_*sl2 )/( 2.*s_ ), // definition from eq. (A.4) in [1]
t1_min = ( w31_*d3+( d3-w31_ )*( d3*w1_-w31_*w2_ )/s_ )/t1_max; // definition from eq. (A.5) in [1]
// FIXME dropped in CDF version
const Kinematics::Limits q2_limits = cuts_.initial_cuts[Cuts::q2];
if ( t1_max > -q2_limits.min() ) { InWarning( Form( "t1max = %f > -q2min = %f", t1_max, -q2_limits.min() ) ); return false; }
if ( t1_min < -q2_limits.max() && q2_limits.hasMax() ) { Debugging( Form( "t1min = %f < -q2max = %f", t1_min, -q2_limits.max() ) ); return false; }
if ( t1_max < -q2_limits.max() && q2_limits.hasMax() ) t1_max = -q2_limits.max();
if ( t1_min > -q2_limits.min() ) t1_min = -q2_limits.min();
/////
// t1, the first photon propagator, is defined here
t1_ = 0.;
double dt1 = 0.;
Map( x(0), t1_min, t1_max, t1_, dt1, "t1" );
// changes wrt mapt1 : dx->-dx
dt1 = -dt1;
DebuggingInsideLoop( Form( "Definition of t1 = %f according to\n\t"
"(t1min, t1max) = (%f, %f)", t1_, t1_min, t1_max ) );
dd4_ = w4_-t1_;
const double d8 = t1_-w2_,
t13 = t1_-w1_-MX2_;
sa1_ = -pow( t1_-w31_, 2 )/4.+w1_*t1_;
if ( sa1_ >= 0. ) { InWarning( Form( "sa1_ = %f >= 0", sa1_ ) ); return false; }
const double sl3 = sqrt( -sa1_ );
// one computes splus and (s2x=s2max)
double splus, s2max;
if ( w1_ != 0. ) {
const double sb =( s_*( t1_-w31_ )+w12_*t13 )/( 2.*w1_ )+MX2_,
sd = sl1_*sl3/w1_,
se =( s_*( t1_*( s_+t13-w2_ )-w2_*w31_ )+MX2_*( w12_*d8+w2_*MX2_ ) )/w1_;
if ( fabs( ( sb-sd )/sd )>=1. ) { splus = sb-sd; s2max = se/splus; }
else { s2max = sb+sd; splus = se/s2max; }
}
else { // 3
s2max = ( s_*(t1_*(s_+d8-MX2_)-w2_*MX2_ )+w2_*MX2_*( w2_+MX2_-t1_ ) )/( ss*t13 );
splus = sig2;
}
// 4
double s2x = s2max;
DebuggingInsideLoop( Form( "s2x = s2max = %f", s2x ) );
if ( n_opt_ < 0 ) { // 5
if ( splus > sig2 ) {
sig2 = splus;
DebuggingInsideLoop( Form( "sig2 truncated to splus = %f", splus ) );
}
if ( n_opt_ < -1 ) { Map( x(2), sig2, s2max, s2_, ds2, "s2" ); }
else { Mapla( t1_, w2_, x(2), sig2, s2max, s2_, ds2 ); } // n_opt_==-1
s2x = s2_;
}
else if ( n_opt_ == 0 ) { s2x = s2_; } // 6
DebuggingInsideLoop( Form( "s2x = %f", s2x ) );
// 7
const double r1 = s2x-d8,
r2 = s2x-d6;
const double rl4 = ( r1*r1-4.*w2_*s2x )*( r2*r2-4.*MY2_*s2x );
if ( rl4 <= 0. ) {
DebuggingInsideLoop( Form( "rl4 = %f <= 0", rl4 ) );
return false;
}
const double sl4 = sqrt( rl4 );
// t2max, t2min definitions from eq. (A.12) and (A.13) in [1]
const double t2_max = w2_+MY2_-( r1*r2+sl4 )/s2x * 0.5,
t2_min = ( w52_*dd4_+( dd4_-w52_ )*( dd4_*w2_-w52_*t1_ )/s2x )/t2_max;
// t2, the second photon propagator, is defined here
t2_ = 0.;
double dt2 = 0.;
Map( x(1), t2_min, t2_max, t2_, dt2, "t2" );
// changes wrt mapt2 : dx->-dx
dt2 = -dt2;
// \f$\delta_6=m_4^2-m_5^2\f$ as defined in Vermaseren's paper
const double tau = t1_-t2_,
r3 = dd4_-t2_,
r4 = w52_-t2_;
DebuggingInsideLoop( Form( "r1 = %f\n\tr2 = %f\n\tr3 = %f\n\tr4 = %f", r1, r2, r3, r4 ) );
const double b = r3*r4-2.*( t1_+w2_ )*t2_,
c = t2_*d6*d8+( d6-d8 )*( d6*w2_-d8*MY2_ );
const double t25 = t2_-w2_-MY2_;
sa2_ = -r4*r4/4.+w2_*t2_;
if ( sa2_ >= 0. ) { InWarning( Form( "sa2_ = %f >= 0", sa2_ ) ); return false; }
const double sl6 = 2.*sqrt( -sa2_ );
g4_ = -r3*r3/4.+t1_*t2_;
if ( g4_ >= 0. ) { InWarning( Form( "g4_ = %f >= 0", g4_ ) ); return false; }
const double sl7 = 2.*sqrt( -g4_ ),
sl5 = sl6*sl7;
double s2p, s2min;
if ( fabs( ( sl5-b )/sl5 ) >= 1. ) {
s2p = ( sl5-b )/t2_ * 0.5;
s2min = c/( t2_*s2p );
}
else { // 8
s2min = ( -sl5-b )/t2_ * 0.5;
s2p = c/( t2_*s2min );
}
// 9
if ( n_opt_ > 1 ) Map( x( 2 ), s2min, s2max, s2_, ds2, "s2" );
else if ( n_opt_ == 1 ) Mapla( t1_, w2_, x( 2 ), s2min, s2max, s2_, ds2 );
const double ap = -0.25*pow( s2_+d8, 2 )+s2_*t1_;
DebuggingInsideLoop( Form( "s2 = %f, s2max = %f, splus = %f", s2_, s2max, splus ) );
if ( w1_!=0. ) dd1_ = -0.25 * ( s2_-s2max ) * ( s2_-splus ) * w1_; // 10
else dd1_ = 0.25 * ( s2_-s2max ) * ss * t13;
// 11
dd2_ = -t2_*( s2_-s2p )*( s2_-s2min ) * 0.25;
DebuggingInsideLoop( Form( "t2 =%f\n\ts2 =%f\n\ts2p=%f\n\ts2min=%f\n\tdd2=%f", t2_, s2_, s2p, s2min, dd2_ ) );
const double yy4 = cos( M_PI*x( 3 ) );
const double dd = dd1_*dd2_;
p12_ = ( s_-w1_-w2_ )*0.5;
const double st = s2_-t1_-w2_;
const double delb = ( 2.*w2_*r3+r4*st )*( 4.*p12_*t1_-( t1_-w31_ )*st )/( 16.*ap );
if ( dd <= 0. ) {
DebuggingInsideLoop( Form( "dd = %e <= 0\n\tdd1 = %e\tdd2 = %e", dd, dd1_, dd2_ ) );
return false;
}
delta_ = delb - yy4*st*sqrt( dd )/ ap * 0.5;
s1_ = t2_+w1_+( 2.*p12_*r3-4.*delta_ )/st;
if ( ap >= 0. ) {
DebuggingInsideLoop( Form( "ap = %f >= 0", ap ) );
return false;
}
jacobian_ = ds2
* dt1
* dt2
* M_PI*M_PI/( 8.*sl1_*sqrt( -ap ) );
DebuggingInsideLoop( Form( "Jacobian = %e", jacobian_ ) );
gram_ = ( 1.-yy4*yy4 )*dd/ap;
p13_ = -t13 * 0.5;
p14_ = (tau+s1_-MX2_) * 0.5;
p25_ = -t25 * 0.5;
/*const double //p15 = (s_+t2_-s1_-w2_)/2.,
//p23 = (s_+t1_-s2_-w1_)/2.,
//p24 = (s2_-tau-MY2_)/2.,
//p34 = (s1_-MX2_-w4_)/2.,
//p35 = (s_+w4_-s1_-s2_)/2.,
//p45 = (s2_-w4_-MY2_)/2.;*/
p1k2_ = ( s1_-t2_-w1_ ) * 0.5;
p2k1_ = st * 0.5;
double s1p, s1m;
if ( w2_ != 0. ) {
const double sbb = ( s_*( t2_-w52_ )-w12_*t25 )/w2_ * 0.5 + MY2_,
sdd = sl1_*sl6/w2_ * 0.5,
see = ( s_*( t2_*( s_+t25-w1_ )-w1_*w52_ )+MY2_*( w1_*MY2_-w12_*( t2_-w1_ ) ) )/w2_;
if ( sbb/sdd >= 0. ) { s1p = sbb+sdd; s1m = see/s1p; }
else { s1m = sbb-sdd; s1p = see/s1m; } // 12
dd3_ = -w2_*( s1p-s1_ )*( s1m-s1_ ) * 0.25; // 13
}
else { // 14
s1p = ( s_*( t2_*( s_-MY2_+t2_-w1_ )-w1_*MY2_ )+w1_*MY2_*( w1_+MY2_-t2_ ) )/( t25*( s_-w12_ ) );
dd3_ = -t25*( s_-w12_ )*( s1p-s1_ ) * 0.25;
}
// 15
//const double acc3 = (s1p-s1_)/(s1p+s1_);
const double ssb = t2_+w1_-r3*( w31_-t1_ )/t1_ * 0.5,
ssd = sl3*sl7/t1_,
sse = ( t2_-w1_ )*( w4_-MX2_ )+( t2_-w4_+w31_ )*( ( t2_-w1_)*MX2_-(w4_-MX2_)*w1_)/t1_;
double s1pp, s1pm;
if (ssb/ssd>=0.) { s1pp = ssb+ssd; s1pm = sse/s1pp; }
else { s1pm = ssb-ssd; s1pp = sse/s1pm; } // 16
// 17
dd4_ = -t1_*( s1_-s1pp )*( s1_-s1pm ) * 0.25;
//const double acc4 = ( s1_-s1pm )/( s1_+s1pm );
dd5_ = dd1_+dd3_+( ( p12_*( t1_-w31_ )*0.5-w1_*p2k1_ )*( p2k1_*( t2_-w52_ )-w2_*r3 )-delta_*( 2.*p12_*p2k1_-w2_*( t1_-w31_ ) ) ) / p2k1_;
return true;
}
bool
GamGamLL::orient()
{
if ( !pickin() or jacobian_ == 0. ) { DebuggingInsideLoop( Form( "Pickin failed! Jacobian = %f", jacobian_ ) ); return false; }
const double re = 0.5 / sqs_;
ep1_ = re*( s_+w12_ );
ep2_ = re*( s_-w12_ );
DebuggingInsideLoop( Form( " re = %e\n\tw12_ = %e", re, w12_ ) );
DebuggingInsideLoop( Form( "Incoming particles' energy = %f, %f", ep1_, ep2_ ) );
p_cm_ = re*sl1_;
de3_ = re*(s2_-MX2_+w12_);
de5_ = re*(s1_-MY2_-w12_);
// Final state energies
const double ep3 = ep1_-de3_,
ep5 = ep2_-de5_;
ec4_ = de3_+de5_;
if ( ec4_ < mc4_ ) { InWarning( Form( "ec4_ = %f < mc4_ = %f\n\t==> de3 = %f, de5 = %f", ec4_, mc4_, de3_, de5_ ) ); return false; }
// What if the protons' momenta are not along the z-axis?
pc4_ = sqrt( ec4_*ec4_-mc4_*mc4_ );
if ( pc4_==0. ) {
InWarning( "pzc4==0" );
return false;
}
const double pp3 = sqrt( ep3*ep3-MX2_ ), pt3 = sqrt( dd1_/s_ )/p_cm_,
pp5 = sqrt( ep5*ep5-MY2_ ), pt5 = sqrt( dd3_/s_ )/p_cm_;
DebuggingInsideLoop( Form( "Central system's energy: E4 = %f\n\t"
" momentum: p4 = %f\n\t"
" invariant mass: m4 = %f\n\t"
"Outgoing particles' energy: E3 = %f\n\t"
" E5 = %f",
ec4_, pc4_, mc4_, ep3, ep5 ) );
const double sin_theta3 = pt3/pp3,
sin_theta5 = pt5/pp5;
DebuggingInsideLoop( Form( "sin_theta3 = %e\n\tsin_theta5 = %e", sin_theta3, sin_theta5 ) );
if (sin_theta3>1.) { InWarning( Form( "sin_theta3 = %e > 1", sin_theta3 ) ); return false; }
if (sin_theta5>1.) { InWarning( Form( "sin_theta5 = %e > 1", sin_theta5 ) ); return false; }
double ct3 = sqrt( 1.-sin_theta3*sin_theta3 ),
ct5 = sqrt( 1.-sin_theta5*sin_theta5 );
if ( ep1_*ep3 < p13_ ) ct3 *= -1.;
if ( ep2_*ep5 > p25_ ) ct5 *= -1.;
DebuggingInsideLoop( Form( "ct3 = %e\n\tct5 = %e", ct3, ct5 ) );
if ( dd5_ < 0. ) { InWarning( Form( "dd5 = %f < 0", dd5_ ) ); return false; }
// Centre of mass system kinematics (theta4 and phi4)
pt4_ = sqrt( dd5_/s_ )/p_cm_;
sin_theta4_ = pt4_/pc4_;
if ( sin_theta4_ > 1. ) { InWarning( Form( "st4 = %f > 1", sin_theta4_ ) ); return false; }
cos_theta4_ = sqrt( 1.-sin_theta4_*sin_theta4_ );
if ( ep1_*ec4_ < p14_ ) cos_theta4_ *= -1.;
al4_ = 1.-cos_theta4_;
be4_ = 1.+cos_theta4_;
if (cos_theta4_<0.) be4_ = sin_theta4_*sin_theta4_/al4_;
else al4_ = sin_theta4_*sin_theta4_/be4_;
DebuggingInsideLoop( Form( "ct4 = %f\n\tal4 = %f, be4 = %f", cos_theta4_, al4_, be4_ ) );
const double rr = sqrt( -gram_/s_ )/( p_cm_*pt4_ );
const double sin_phi3 = rr / pt3,
sin_phi5 = -rr / pt5;
if ( fabs( sin_phi3 ) > 1. ) { InWarning( Form( "sin(phi_3) = %e while it must be in [-1 ; 1]", sin_phi3 ) ); return false; }
if ( fabs( sin_phi5 ) > 1. ) { InWarning( Form( "sin(phi_5) = %e while it must be in [-1 ; 1]", sin_phi5 ) ); return false; }
const double cos_phi3 = -sqrt( 1.-sin_phi3*sin_phi3 ),
cos_phi5 = -sqrt( 1.-sin_phi5*sin_phi5 );
p3_lab_ = Particle::Momentum( pp3*sin_theta3*cos_phi3, pp3*sin_theta3*sin_phi3, pp3*ct3, ep3 );
p5_lab_ = Particle::Momentum( pp5*sin_theta5*cos_phi5, pp5*sin_theta5*sin_phi5, pp5*ct5, ep5 );
const double a1 = p3_lab_.px()-p5_lab_.px();
DebuggingInsideLoop( Form( "Kinematic quantities\n\t"
"cos(theta3) = %1.4f\tsin(theta3) = %1.4f\tcos( phi3 ) = %1.4f\tsin( phi3 ) = %1.4f\n\t"
"cos(theta4) = %1.4f\tsin(theta4) = %1.4f\n\t"
"cos(theta5) = %1.4f\tsin(theta5) = %1.4f\tcos( phi5 ) = %1.4f\tsin( phi5 ) = %1.4f\n\t"
"a1 = %f",
ct3, sin_theta3, cos_phi3, sin_phi3,
cos_theta4_, sin_theta4_,
ct5, sin_theta5, cos_phi5, sin_phi5, a1 ) );
if ( fabs( pt4_+p3_lab_.px()+p5_lab_.px() ) < fabs( fabs( a1 )-pt4_ ) ) {
DebuggingInsideLoop( Form( "|pp4+pp3*cos(phi3)+pp5*cos(phi5)| < | |a1|-pp4 |\n\t"
"pp4 = %f\tpp5 = %f\n\t"
"cos(phi3) = %f\tcos(phi5) = %f"
"a1 = %f",
pt4_, pt5, cos_phi3, cos_phi5, a1 ) );
return true;
}
if ( a1 < 0. ) p5_lab_.setP( 0, -p5_lab_.px() );
else p3_lab_.setP( 0, -p3_lab_.px() );
return true;
}
double
GamGamLL::computeOutgoingPrimaryParticlesMasses( double x, double outmass, double lepmass, double& dw )
{
const double mx0 = Particle::massFromPDGId( Particle::Proton )+Particle::massFromPDGId( Particle::PiPlus ); // 1.07
const Kinematics::Limits mx_limits = cuts_.remnant_cuts[Cuts::mass];
const double wx2min = pow( std::max( mx0, mx_limits.min() ), 2 ),
wx2max = pow( std::min( sqs_-outmass-2.*lepmass, mx_limits.max() ), 2 );
double mx2 = 0., dmx2 = 0.;
Map( x, wx2min, wx2max, mx2, dmx2, "mx2" );
DebuggingInsideLoop( Form( "mX^2 in range (%f, %f), x = %f\n\t"
"mX^2 = %f, d(mX^2) = %f\n\t"
"mX = %f, d(mX) = %f", wx2min, wx2max, x, mx2, dmx2, sqrt( mx2 ), sqrt( dmx2 ) ) );
dw = sqrt( dmx2 );
return sqrt( mx2 );
}
void
GamGamLL::beforeComputeWeight()
{
if ( !GenericProcess::is_point_set_ ) return;
const Particle& p1 = event_->getOneByRole( Particle::IncomingBeam1 ),
&p2 = event_->getOneByRole( Particle::IncomingBeam2 );
ep1_ = p1.energy();
ep2_ = p2.energy();
Ml2_ = event_->getByRole( Particle::CentralSystem )[0].mass2();
switch ( cuts_.mode ) {
case Kinematics::ElectronProton: default:
{ InError( "Case not yet supported!" ); } break;
case Kinematics::ElasticElastic:
dw31_ = dw52_ = 0.; break;
case Kinematics::InelasticElastic: {
const double m = computeOutgoingPrimaryParticlesMasses( x( 7 ), p1.mass(), sqrt( Ml2_ ), dw31_ );
event_->getOneByRole( Particle::OutgoingBeam1 ).setMass( m );
event_->getOneByRole( Particle::OutgoingBeam2 ).setMass( Particle::massFromPDGId( p2.pdgId() ) );
} break;
case Kinematics::ElasticInelastic: {
const double m = computeOutgoingPrimaryParticlesMasses( x( 7 ), p2.mass(), sqrt( Ml2_ ), dw52_ );
event_->getOneByRole( Particle::OutgoingBeam1 ).setMass( Particle::massFromPDGId( p1.pdgId() ) );
event_->getOneByRole( Particle::OutgoingBeam2 ).setMass( m );
} break;
case Kinematics::InelasticInelastic: {
const double mx = computeOutgoingPrimaryParticlesMasses( x( 7 ), p2.mass(), sqrt( Ml2_ ), dw31_ );
event_->getOneByRole( Particle::OutgoingBeam1 ).setMass( mx );
const double my = computeOutgoingPrimaryParticlesMasses( x( 8 ), p1.mass(), sqrt( Ml2_ ), dw52_ );
event_->getOneByRole( Particle::OutgoingBeam2 ).setMass( my );
} break;
}
MX_ = event_->getOneByRole( Particle::OutgoingBeam1 ).mass();
MY_ = event_->getOneByRole( Particle::OutgoingBeam2 ).mass();
MX2_ = MX_*MX_;
MY2_ = MY_*MY_;
}
double
GamGamLL::computeWeight()
{
if ( !is_outgoing_state_set_ ) { InWarning( "Output state not set!" ); return 0.; }
DebuggingInsideLoop( Form( "sqrt(s)=%f\n\tm(X1)=%f\tm(X2)=%f", sqs_, MX_, MY_ ) );
Kinematics::Limits& w_limits = cuts_.initial_cuts[Cuts::w];
if ( !w_limits.hasMax() ) w_limits.max() = s_;
// The minimal energy for the central system is its outgoing leptons' mass energy (or wmin_ if specified)
const double wmin = std::max( 4.*Ml2_, w_limits.min() );
// The maximal energy for the central system is its CM energy with the outgoing particles' mass energy substracted (or _wmax if specified)
const double wmax = std::min( pow( sqs_-MX_-MY_, 2 ), w_limits.max() );
DebuggingInsideLoop( Form( "wmin = %f\n\twmax = %f\n\twmax/wmin = %f", wmin, wmax, wmax/wmin ) );
// compute the two-photon energy for this point
w4_ = 0.;
double dw4 = 0.;
Map( x( 4 ), wmin, wmax, w4_, dw4, "w4" );
Map( x( 4 ), wmin, wmax, w4_, dw4, "w4" );
mc4_ = sqrt( w4_ );
DebuggingInsideLoop( Form( "Computed value for w4 = %f -> mc4 = %f", w4_, mc4_ ) );
if ( !orient() ) return 0.;
if ( jacobian_ == 0. ) { InWarning( Form( "dj = %f", jacobian_ ) ); return 0.; }
if ( t1_>0. ) { InWarning( Form( "t1 = %f > 0", t1_ ) ); return 0.; }
if ( t2_>0. ) { InWarning( Form( "t2 = %f > 0", t2_ ) ); return 0.; }
const double ecm6 = w4_ / ( 2.*mc4_ ),
pp6cm = sqrt( ecm6*ecm6-Ml2_ );
jacobian_ *= dw4*pp6cm/( mc4_*Constants::sconstb*s_ );
// Let the most obscure part of this code begin...
const double e1mp1 = w1_ / ( ep1_+p_cm_ ),
e3mp3 = MX2_ / ( p3_lab_.energy()+p3_lab_.p() );
const double al3 = pow( sin( p3_lab_.theta() ), 2 )/( 1.+( p3_lab_.theta() ) );
// 2-photon system kinematics ?!
const double eg = ( w4_+t1_-t2_ )/( 2.*mc4_ );
double pg = sqrt( eg*eg-t1_ );
const double pgx = -p3_lab_.px()*cos_theta4_-sin_theta4_*( de3_-e1mp1 + e3mp3 + p3_lab_.p()*al3 ),
pgy = -p3_lab_.py(),
pgz = mc4_*de3_/( ec4_+pc4_ )-ec4_*de3_*al4_/mc4_-p3_lab_.px()*ec4_*sin_theta4_/mc4_+ec4_*cos_theta4_/mc4_*( p3_lab_.p()*al3+e3mp3-e1mp1 );
DebuggingInsideLoop( Form( "pg3 = (%f, %f, %f)\n\t"
"pg3^2 = %f",
pgx, pgy, pgz,
sqrt( pgx*pgx+pgy*pgy+pgz*pgz ) ) );
const double pgp = sqrt( pgx*pgx + pgy*pgy ), // outgoing proton (3)'s transverse momentum
pgg = sqrt( pgp*pgp + pgz*pgz ); // outgoing proton (3)'s momentum
if ( pgg > pgp*0.9 && pgg > pg ) { pg = pgg; } //FIXME ???
// Phi angle for the 2-photon system ?!
const double cpg = pgx/pgp,
spg = pgy/pgp;
// Theta angle for the 2-photon system ?!
const double stg = pgp/pg;
const int theta_sign = ( pgz>0. ) ? 1 : -1;
const double ctg = theta_sign*sqrt( 1.-stg*stg );
double xx6 = x( 5 );
const double amap = 0.5 * (w4_-t1_-t2_),
bmap = 0.5 * sqrt( ( pow( w4_-t1_-t2_, 2 )-4.*t1_*t2_ )*( 1.-4.*Ml2_/w4_ ) ),
ymap = ( amap+bmap )/( amap-bmap ),
beta = pow( ymap, 2.*xx6-1. );
xx6 = 0.5 * ( 1. + amap/bmap*( beta-1. )/( beta+1. ) );
xx6 = std::max( 0., std::min( xx6, 1. ) ); // xx6 in [0., 1.]
DebuggingInsideLoop( Form( "amap = %f\n\tbmap = %f\n\tymap = %f\n\tbeta = %f", amap, bmap, ymap, beta ) );
// 3D rotation of the first outgoing lepton wrt the CM system
const double theta6cm = acos( 1.-2.*xx6 );
// match the Jacobian
jacobian_ *= ( amap+bmap*cos( theta6cm ) );
jacobian_ *= ( amap-bmap*cos( theta6cm ) );
jacobian_ /= amap;
jacobian_ /= bmap;
jacobian_ *= log( ymap );
jacobian_ *= 0.5;
DebuggingInsideLoop( Form( "Jacobian = %e", jacobian_ ) );
DebuggingInsideLoop( Form( "ctcm6 = %f\n\tstcm6 = %f", cos( theta6cm ), sin( theta6cm ) ) );
const double phi6cm = 2.*M_PI*x( 6 );
// First outgoing lepton's 3-momentum in the centre of mass system
Particle::Momentum p6cm = Particle::Momentum::fromPThetaPhi( pp6cm, theta6cm, phi6cm );
DebuggingInsideLoop( Form( "p3cm6 = (%f, %f, %f)", p6cm.px(), p6cm.py(), p6cm.pz() ) );
const double h1 = stg*p6cm.pz()+ctg*p6cm.px();
const double pc6z = ctg*p6cm.pz()-stg*p6cm.px(), pc6x = cpg*h1-spg*p6cm.py();
const double qcx = 2.*pc6x, qcz = 2.*pc6z;
// qcy == QCY is never defined
const double el6 = ( ec4_*ecm6+pc4_*pc6z ) / mc4_,
h2 = ( ec4_*pc6z+pc4_*ecm6 ) / mc4_;
DebuggingInsideLoop( Form( "h1 = %f\n\th2 = %f", h1, h2 ) );
// first outgoing lepton's 3-momentum
const double p6x = cos_theta4_*pc6x+sin_theta4_*h2,
p6y = cpg*p6cm.py()+spg*h1,
p6z = cos_theta4_*h2-sin_theta4_*pc6x;
// first outgoing lepton's kinematics
p6_cm_ = Particle::Momentum( p6x, p6y, p6z, el6 );
DebuggingInsideLoop( Form( "E6(cm) = %f\n\tP6(cm) = (%f, %f, %f)", el6, p6x, p6y, p6z ) );
const double hq = ec4_*qcz/mc4_;
const Particle::Momentum qve(
cos_theta4_*qcx+sin_theta4_*hq,
2.*p6y,
cos_theta4_*hq-sin_theta4_*qcx,
pc4_*qcz/mc4_ // energy
);
// Available energy for the second lepton is the 2-photon system's energy with the first lepton's energy removed
const double el7 = ec4_-el6;
DebuggingInsideLoop( Form( "Outgoing kinematics\n\t"
" first outgoing lepton: p = %f, E = %f\n\t"
"second outgoing lepton: p = %f, E = %f",
p6_cm_.p(), p6_cm_.energy(), p7_cm_.p(), p6_cm_.energy() ) );
// Second outgoing lepton's 3-momentum
const double p7x = -p6x + pt4_,
p7y = -p6y,
p7z = -p6z + pc4_*cos_theta4_;
// second outgoing lepton's kinematics
p7_cm_ = Particle::Momentum( p7x, p7y, p7z, el7 );
//p6_cm_ = Particle::Momentum(pl6*st6*cp6, pl6*st6*sp6, pl6*ct6, el6);
//p7_cm_ = Particle::Momentum(pl7*st7*cp7, pl7*st7*sp7, pl7*ct7, el7);
q1dq_ = eg*( 2.*ecm6-mc4_ )-2.*pg*p6cm.pz();
q1dq2_ = ( w4_-t1_-t2_ ) * 0.5;
const double phi3 = p3_lab_.phi(), cos_phi3 = cos( phi3 ), sin_phi3 = sin( phi3 ),
phi5 = p5_lab_.phi(), cos_phi5 = cos( phi5 ), sin_phi5 = sin( phi5 );
bb_ = t1_*t2_+( w4_*pow( sin( theta6cm ), 2 ) + 4.*Ml2_*pow( cos( theta6cm ), 2 ) )*pg*pg;
const double c1 = p3_lab_.pt() * ( qve.px()*sin_phi3 - qve.py()*cos_phi3 ),
c2 = p3_lab_.pt() * ( qve.pz()*ep1_ - qve.energy() *p_cm_ ),
c3 = ( w31_*ep1_*ep1_ + 2.*w1_*de3_*ep1_ - w1_*de3_*de3_ + p3_lab_.pt2()*ep1_*ep1_ ) / ( p3_lab_.energy()*p_cm_ + p3_lab_.pz()*ep1_ );
const double b1 = p5_lab_.pt() * ( qve.px()*sin_phi5 - qve.py()*cos_phi5 ),
b2 = p5_lab_.pt() * ( qve.pz()*ep2_ + qve.energy() *p_cm_ ),
b3 = ( w52_*ep2_*ep2_ + 2.*w2_*de5_*ep2_ - w2_*de5_*de5_ + p5_lab_.pt2()*ep2_*ep2_ ) / ( ep2_*p5_lab_.pz() - p5_lab_.energy()*p_cm_ );
const double r12 = c2*sin_phi3 + qve.py()*c3,
r13 = -c2*cos_phi3 - qve.px()*c3;
const double r22 = b2*sin_phi5 + qve.py()*b3,
r23 = -b2*cos_phi5 - qve.px()*b3;
epsi_ = p12_*c1*b1 + r12*r22 + r13*r23;
g5_ = w1_*c1*c1 + r12*r12 + r13*r13;
g6_ = w2_*b1*b1 + r22*r22 + r23*r23;
a5_ = -( qve.px()*cos_phi3 + qve.py()*sin_phi3 )*p3_lab_.pt()*p1k2_-( ep1_*qve.energy()-p_cm_*qve.pz() )*( cos_phi3*cos_phi5 + sin_phi3*sin_phi5 )*p3_lab_.pt()*p5_lab_.pt() + ( de5_*qve.pz()+qve.energy()*( p_cm_+p5_lab_.pz() ) )*c3;
a6_ = -( qve.px()*cos_phi5 + qve.py()*sin_phi5 )*p5_lab_.pt()*p2k1_-( ep2_*qve.energy()+p_cm_*qve.pz() )*( cos_phi3*cos_phi5 + sin_phi3*sin_phi5 )*p3_lab_.pt()*p5_lab_.pt() + ( de3_*qve.pz()-qve.energy()*( p_cm_-p3_lab_.pz() ) )*b3;
DebuggingInsideLoop( Form( "a5 = %f\n\ta6 = %f", a5_, a6_ ) );
////////////////////////////////////////////////////////////////
// END of GAMGAMLL subroutine in the FORTRAN version
////////////////////////////////////////////////////////////////
const Particle::Momentum cm = event_->getOneByRole( Particle::IncomingBeam1 ).momentum() + event_->getOneByRole( Particle::IncomingBeam2 ).momentum();
////////////////////////////////////////////////////////////////
// INFO from f.f
////////////////////////////////////////////////////////////////
const double gamma = cm.energy() / sqs_, betgam = cm.pz() / sqs_;
//--- kinematics computation for both leptons
p6_cm_.betaGammaBoost( gamma, betgam );
p7_cm_.betaGammaBoost( gamma, betgam );
//--- cut on mass of final hadronic system (MX/Y)
const Kinematics::Limits mx_limits = cuts_.remnant_cuts[Cuts::mass];
if ( cuts_.mode == Kinematics::InelasticElastic || cuts_.mode == Kinematics::InelasticInelastic ) {
if ( mx_limits.hasMin() && MX_ < mx_limits.min() ) return 0.;
if ( mx_limits.hasMax() && MX_ > mx_limits.max() ) return 0.;
}
if ( cuts_.mode == Kinematics::ElasticInelastic || cuts_.mode == Kinematics::InelasticInelastic ) {
if ( mx_limits.hasMin() && MY_ < mx_limits.min() ) return 0.;
if ( mx_limits.hasMax() && MY_ > mx_limits.max() ) return 0.;
}
//--- cut on the proton's Q2 (first photon propagator T1)
const Kinematics::Limits q2_limits = cuts_.initial_cuts[Cuts::q2];
if ( q2_limits.hasMax() && t1_<-q2_limits.max() ) return 0;
if ( q2_limits.hasMin() && t1_>-q2_limits.min() ) return 0.;
//--- cuts on outgoing leptons' kinematics
const Kinematics::Limits m_limits = cuts_.central_cuts[Cuts::mass_sum];
if ( m_limits.hasMin() && ( p6_cm_+p7_cm_ ).mass() < m_limits.min() ) return 0.;
if ( m_limits.hasMax() && ( p6_cm_+p7_cm_ ).mass() > m_limits.max() ) return 0.;
//----- cuts on the individual leptons
const Kinematics::Limits pt_limits = cuts_.central_cuts[Cuts::pt_single];
if ( pt_limits.hasMin() && ( p6_cm_.pt() < pt_limits.min() || p7_cm_.pt() < pt_limits.min() ) ) return 0.;
if ( pt_limits.hasMax() && ( p6_cm_.pt() > pt_limits.max() || p7_cm_.pt() > pt_limits.max() ) ) return 0.;
const Kinematics::Limits energy_limits = cuts_.central_cuts[Cuts::energy_single];
if ( energy_limits.hasMin() && ( p6_cm_.energy() < energy_limits.min() || p7_cm_.energy() < energy_limits.min() ) ) return 0.;
if ( energy_limits.hasMax() && ( p6_cm_.energy() > energy_limits.max() || p7_cm_.energy() > energy_limits.max() ) ) return 0.;
const Kinematics::Limits eta_limits = cuts_.central_cuts[Cuts::eta_single];
if ( eta_limits.hasMin() && ( p6_cm_.eta() < eta_limits.min() || p7_cm_.eta() < eta_limits.min() ) ) return 0.;
if ( eta_limits.hasMax() && ( p6_cm_.eta() > eta_limits.max() || p7_cm_.eta() > eta_limits.max() ) ) return 0.;
//--- compute the structure functions factors
switch ( cuts_.mode ) { // inherited from CDF version
case Kinematics::ElectronProton: default: jacobian_ *= periPP( 1, 2 ); break; // ep case
case Kinematics::ElasticElastic: jacobian_ *= periPP( 2, 2 ); break; // elastic case
case Kinematics::InelasticElastic: jacobian_ *= periPP( 3, 2 )*( dw31_*dw31_ ); break;
case Kinematics::ElasticInelastic: jacobian_ *= periPP( 3, 2 )*( dw52_*dw52_ ); break; // single-dissociative case
case Kinematics::InelasticInelastic: jacobian_ *= periPP( 3, 3 )*( dw31_*dw31_ )*( dw52_*dw52_ ); break; // double-dissociative case
}
//--- compute the event weight using the Jacobian
return Constants::GeV2toBarn*jacobian_;
}
void
GamGamLL::fillKinematics( bool )
{
const Particle::Momentum cm = event_->getOneByRole( Particle::IncomingBeam1 ).momentum() + event_->getOneByRole( Particle::IncomingBeam2 ).momentum();
const double gamma = cm.energy()/sqs_, betgam = cm.pz()/sqs_;
Particle::Momentum plab_ip1 = Particle::Momentum( 0., 0., p_cm_, ep1_ ).betaGammaBoost( gamma, betgam );
Particle::Momentum plab_ip2 = Particle::Momentum( 0., 0., -p_cm_, ep2_ ).betaGammaBoost( gamma, betgam );
p3_lab_.betaGammaBoost( gamma, betgam );
p5_lab_.betaGammaBoost( gamma, betgam );
//----- needed to parametrise a random rotation around z-axis
const int rany = ( rand() >= RAND_MAX*0.5 ) ? 1 : -1,
ransign = ( rand() >= RAND_MAX*0.5 ) ? 1 : -1;
const double ranphi = ( rand()/RAND_MAX*0.5 )*2.*M_PI;
Particle::Momentum plab_ph1 = ( plab_ip1-p3_lab_ ).rotatePhi( ranphi, rany );
Particle::Momentum plab_ph2 = ( plab_ip2-p5_lab_ ).rotatePhi( ranphi, rany );
p3_lab_.rotatePhi( ranphi, rany );
p5_lab_.rotatePhi( ranphi, rany );
p6_cm_.rotatePhi( ranphi, rany );
p7_cm_.rotatePhi( ranphi, rany );
/*if ( symmetrise_ && rand() >= .5*RAND_MAX ) {
p6_cm_.mirrorZ();
p7_cm_.mirrorZ();
}*/
//----- incoming protons
event_->getOneByRole( Particle::IncomingBeam1 ).setMomentum( plab_ip1 );
event_->getOneByRole( Particle::IncomingBeam2 ).setMomentum( plab_ip2 );
//----- first outgoing proton
Particle& op1 = event_->getOneByRole( Particle::OutgoingBeam1 );
op1.setMomentum( p3_lab_ );
switch ( cuts_.mode ) {
case Kinematics::ElasticElastic:
case Kinematics::ElasticInelastic:
default:
op1.setStatus( Particle::FinalState ); // stable proton
break;
case Kinematics::InelasticElastic:
case Kinematics::InelasticInelastic:
op1.setStatus( Particle::Undecayed ); // fragmenting remnants
op1.setMass( MX_ );
break;
}
//----- second outgoing proton
Particle& op2 = event_->getOneByRole( Particle::OutgoingBeam2 );
op2.setMomentum( p5_lab_ );
switch ( cuts_.mode ) {
case Kinematics::ElasticElastic:
case Kinematics::InelasticElastic:
default:
op2.setStatus( Particle::FinalState ); // stable proton
break;
case Kinematics::ElasticInelastic:
case Kinematics::InelasticInelastic:
op2.setStatus( Particle::Undecayed ); // fragmenting remnants
op2.setMass( MY_ );
break;
}
//----- first incoming photon
Particle& ph1 = event_->getOneByRole( Particle::Parton1 );
ph1.setMomentum( plab_ph1 );
ph1.setStatus( Particle::Incoming );
//----- second incoming photon
Particle& ph2 = event_->getOneByRole( Particle::Parton2 );
ph2.setMomentum( plab_ph2 );
ph2.setStatus( Particle::Incoming );
Particles& central_system = event_->getByRole( Particle::CentralSystem );
//----- first outgoing lepton
Particle& ol1 = central_system[0];
ol1.setPdgId( ol1.pdgId(), ransign );
ol1.setMomentum( p6_cm_ );
ol1.setStatus( Particle::FinalState );
//----- second outgoing lepton
Particle& ol2 = central_system[1];
ol2.setPdgId( ol2.pdgId(), -ransign );
ol2.setMomentum( p7_cm_ );
ol2.setStatus( Particle::FinalState );
//----- intermediate two-lepton system
event_->getOneByRole( Particle::Intermediate ).setMomentum( p6_cm_+p7_cm_ );
}
double
GamGamLL::periPP( int nup_, int ndown_ )
{
DebuggingInsideLoop( Form( " Nup = %d\n\tNdown = %d", nup_, ndown_ ) );
FormFactors fp1, fp2;
GenericProcess::formFactors( -t1_, -t2_, fp1, fp2 );
DebuggingInsideLoop( Form( "u1 = %f\n\tu2 = %f\n\tv1 = %f\n\tv2 = %f", fp1.FM, fp1.FE, fp2.FM, fp2.FE ) );
const double qqq = q1dq_*q1dq_,
qdq = 4.*Ml2_-w4_;
const double t11 = 64. *( bb_*( qqq-g4_-qdq*( t1_+t2_+2.*Ml2_ ) )-2.*( t1_+2.*Ml2_ )*( t2_+2.*Ml2_ )*qqq ) * t1_*t2_, // magnetic-magnetic
t12 = 128.*( -bb_*( dd2_+g6_ )-2.*( t1_+2.*Ml2_ )*( sa2_*qqq+a6_*a6_ ) ) * t1_, // electric-magnetic
t21 = 128.*( -bb_*( dd4_+g5_ )-2.*( t2_+2.*Ml2_ )*( sa1_*qqq+a5_*a5_ ) ) * t2_, // magnetic-electric
t22 = 512.*( bb_*( delta_*delta_-gram_ )-pow( epsi_-delta_*( qdq+q1dq2_ ), 2 )-sa1_*a6_*a6_-sa2_*a5_*a5_-sa1_*sa2_*qqq ); // electric-electric
const double peripp = ( fp1.FM*fp2.FM*t11 + fp1.FE*fp2.FM*t21 + fp1.FM*fp2.FE*t12 + fp1.FE*fp2.FE*t22 ) / pow( 2.*t1_*t2_*bb_, 2 );
DebuggingInsideLoop( Form( "t11 = %5.2f\tt12 = %5.2f\n\t"
"t21 = %5.2f\tt22 = %5.2f\n\t"
"=> PeriPP = %e", t11, t12, t21, t22, peripp ) );
return peripp;
}
diff --git a/CepGen/Processes/GamGamLL.h b/CepGen/Processes/GamGamLL.h
index 9082ce6..c5c758a 100644
--- a/CepGen/Processes/GamGamLL.h
+++ b/CepGen/Processes/GamGamLL.h
@@ -1,207 +1,205 @@
#ifndef CepGen_Processes_GamGamLL_h
#define CepGen_Processes_GamGamLL_h
#define PROCESS_DESCRIPTION "Two-photon production of lepton pairs"
#include <iomanip>
#include <algorithm>
#include <string>
#include <cmath>
#include <ctime>
#include <cstdlib>
#include <map>
#include "GenericProcess.h"
#include "CepGen/Physics/FormFactors.h"
namespace CepGen
{
namespace Process
{
/**
* Full class of methods and objects to compute the full analytic matrix element
* \cite Vermaseren1983347 for the \f$\gamma\gamma\to\ell^{+}\ell^{-}\f$ process
* according to a set of kinematic constraints provided for the incoming and
* outgoing particles (the Kinematics object).
* The particle roles in this process are defined as following: \n
* \image latex lpair_kinematics.pdf Detailed particle roles in the two-photon process as defined by the @a GamGamLL object. The incoming protons/electrons are denoted by a role 1, and 2, as the outgoing protons/protons remnants/ electrons carry the indices 3 and 5. The two outgoing leptons have the roles 6 and 7, while the lepton/antilepton distinction is done randomly (thus, the arrow convention is irrelevant here).
*
* The \a f function created by this Process child has its \a _ndim -dimensional
* coordinates mapped as :
* - 0 = \f$t_1\f$, first incoming photon's virtuality
* - 1 = \f$t_2\f$, second incoming photon's virtuality
* - 2 = \f$s_2\f$ mapping
* - 3 = yy4 = \f$\cos\left(\pi x_3\right)\f$ definition
* - 4 = \f$w_4\f$, the two-photon system's invariant mass
* - 5 = xx6 = \f$\frac{1}{2}\left(1-\cos\theta^\text{CM}_6\right)\f$ definition (3D rotation of the first outgoing lepton with respect to the two-photon centre-of-mass system). If the @a nm_ optimisation flag is set this angle coefficient value becomes
* \f[\frac{1}{2}\left(\frac{a_\text{map}}{b_\text{map}}\frac{\beta-1}{\beta+1}+1\right)\f]
* with \f$a_\text{map}=\frac{1}{2}\left(w_4-t_1-t_2\right)\f$, \f$b_\text{map}=\frac{1}{2}\sqrt{\left(\left(w_4-t_1-t_2\right)^2-4t_1t_2\right)\left(1-4\frac{w_6}{w_4}\right)}\f$, and \f$\beta=\left(\frac{a_\text{map}+b_\text{map}}{a_\text{map}-b_\text{map}}\right)^{2x_5-1}\f$
* and the \a fJacobian element is scaled by a factor \f$\frac{1}{2}\frac{\left(a_\text{map}^2-b_\text{map}^2\cos^2\theta^\text{CM}_6\right)}{a_\text{map}b_\text{map}}\log\left(\frac{a_\text{map}+b_\text{map}}{a_\text{map}-b_\text{map}}\right)\f$
* - 6 = _phicm6_, or \f$\phi_6^\text{CM}\f$ the rotation angle of the dilepton system in the centre-of-mass
* system
* - 7 = \f$x_q\f$, \f$w_X\f$ mappings, as used in the single- and double-dissociative
* cases only
* \brief Compute the matrix element for a CE \f$\gamma\gamma\to\ell^{+}\ell^{-}\f$
* process
*/
class GamGamLL : public GenericProcess
{
public:
/// Class constructor ; set the mandatory parameters before integration and events generation
/// \param[in] nopt Optimisation (legacy from LPAIR)
GamGamLL( int nopt=0 );
void addEventContent();
void beforeComputeWeight();
/// Compute the process' weight for the given point
/// \return \f$\mathrm d\sigma(\mathbf x)(\gamma\gamma\to\ell^{+}\ell^{-})\f$,
/// the differential cross-section for the given point in the phase space.
double computeWeight();
unsigned int numDimensions( const Kinematics::ProcessMode& ) const;
void fillKinematics( bool );
/// Compute the ougoing proton remnant mass
/// \param[in] x A random number (between 0 and 1)
/// \param[in] outmass The maximal outgoing particles' invariant mass
/// \param[in] lepmass The outgoing leptons' mass
/// \param[out] dw The size of the integration bin
/// \return Mass of the outgoing proton remnant
double computeOutgoingPrimaryParticlesMasses( double x, double outmass, double lepmass, double& dw );
/// Set all the kinematic variables for the outgoing proton remnants, and prepare the hadronisation
/// \param[in] part_ Particle to "prepare" for the hadronisation to be performed
void prepareHadronisation( Particle *part_ );
private:
/**
* Calculate energies and momenta of the
* 1st, 2nd (resp. the "proton-like" and the "electron-like" incoming particles),
* 3rd (the "proton-like" outgoing particle),
* 4th (the two-photons central system), and
* 5th (the "electron-like" outgoing particle) particles in the overall centre-of-mass frame.
* \brief Energies/momenta computation for the various particles, in the CM system
* \return Success state of the operation
*/
bool orient();
/**
* Compute the expression of the matrix element squared for the \f$\gamma\gamma\rightarrow\ell^{+}\ell^{-}\f$ process.
* It returns the value of the convolution of the form factor or structure functions with the central two-photons matrix element squared.
* \brief Computes the matrix element squared for the requested process
* \return Full matrix element for the two-photon production of a pair of spin\f$-\frac{1}{2}-\f$point particles.
* It is noted as \f[
* M = \frac{1}{4bt_1 t_2}\sum_{i=1}^2\sum_{j=1}^2 u_i v_j t_{ij} = \frac{1}{4}\frac{u_1 v_1 t_{11}+u_2 v_1 t_{21}+u_1 v_2 t_{12}+u_2 v_2 t_{22}}{t_1 t_2 b}
* \f] where \f$b\f$ = \a bb_ is defined in \a ComputeWeight as : \f[
* b = t_1 t_2+\left(w_{\gamma\gamma}\sin^2{\theta^\text{CM}_6}+4m_\ell\cos^2{\theta^\text{CM}_6}\right) p_g^2
* \f]
*/
double periPP( int, int );
/**
* Describe the kinematics of the process \f$p_1+p_2\to p_3+p_4+p_5\f$ in terms of Lorentz-invariant variables.
* These variables (along with others) will then be fed into the \a PeriPP method (thus are essential for the evaluation of the full matrix element).
* \return Success state of the operation
*/
bool pickin();
/// Internal switch for the optimised code version (LPAIR legacy ; unimplemented here)
int n_opt_;
/// squared mass of the first proton-like outgoing particle
double MX2_;
/// squared mass of the second proton-like outgoing particle
double MY2_;
/// squared mass of the outgoing leptons
double Ml2_;
/// energy of the first proton-like incoming particle
double ep1_;
/// energy of the second proton-like incoming particle
double ep2_;
double p_cm_;
/// \f$\delta_2=m_1^2-m_2^2\f$ as defined in Vermaseren's paper
/// \cite Vermaseren1983347 for the full definition of this quantity
double w12_;
/// \f$\delta_1=m_3^2-m_1^2\f$ as defined in Vermaseren's paper
/// \cite Vermaseren1983347 for the full definition of this quantity
double w31_;
double dw31_;
/// \f$\delta_4=m_5^2-m_2^2\f$ as defined in Vermaseren's paper
/// \cite Vermaseren1983347 for the full definition of this quantity
double w52_;
double dw52_;
/// energy of the two-photon central system
double ec4_;
/// 3-momentum norm of the two-photon central system
double pc4_;
/// mass of the two-photon central system
double mc4_;
/// squared mass of the two-photon central system
double w4_;
/// \f$p_{12} = \frac{1}{2}\left(s-m_{p_1}^2-m_{p_2}^2\right)\f$
double p12_;
double p1k2_, p2k1_;
/// \f$p_{13} = -\frac{1}{2}\left(t_1-m_{p_1}^2-m_{p_3}^2\right)\f$
double p13_;
double p14_, p25_;
double q1dq_, q1dq2_;
double s1_, s2_;
double epsi_;
double g5_, g6_;
double a5_, a6_;
double bb_;
double gram_;
double dd1_, dd2_, dd3_;
/// \f$\delta_5=m_4^2-t_1\f$ as defined in Vermaseren's paper
/// \cite Vermaseren1983347 for the full definition of this quantity
double dd4_;
double dd5_;
/**
* Invariant used to tame divergences in the matrix element computation. It is defined as
* \f[\Delta = \left(p_1\cdot p_2\right)\left(q_1\cdot q_2\right)-\left(p_1\cdot q_2\right)\left(p_2\cdot q_1\right)\f]
* with \f$p_i, q_i\f$ the 4-momenta associated to the incoming proton-like particle and to the photon emitted from it.
*/
double delta_;
double g4_;
double sa1_, sa2_;
double sl1_;
/// cosine of the polar angle for the two-photons centre-of-mass system
double cos_theta4_;
/// sine of the polar angle for the two-photons centre-of-mass system
double sin_theta4_;
double al4_;
double be4_;
double de3_, de5_;
double pt4_;
/// Kinematics of the first incoming proton
Particle::Momentum p1_lab_;
/// Kinematics of the second incoming proton
Particle::Momentum p2_lab_;
/// Kinematics of the first outgoing proton
Particle::Momentum p3_lab_;
/// Kinematics of the two-photon system (in the two-proton CM)
Particle::Momentum p4_lab_;
/// Kinematics of the second outgoing proton
Particle::Momentum p5_lab_;
/// Kinematics of the first outgoing lepton (in the two-proton CM)
Particle::Momentum p6_cm_;
/// Kinematics of the second outgoing lepton (in the two-proton CM)
Particle::Momentum p7_cm_;
double jacobian_;
-
- double cot_theta1_, cot_theta2_;
};
}
}
#endif
diff --git a/CepGen/Processes/PPtoLL.cpp b/CepGen/Processes/PPtoLL.cpp
index 4d0d641..32f6dff 100644
--- a/CepGen/Processes/PPtoLL.cpp
+++ b/CepGen/Processes/PPtoLL.cpp
@@ -1,412 +1,412 @@
#include "PPtoLL.h"
#include <assert.h>
using namespace CepGen::Process;
PPtoLL::PPtoLL() :
GenericKTProcess( "pptoll", "gamma,gamma->l+,l-", 4, Particle::Photon, Particle::Muon )
{}
void
PPtoLL::prepareKTKinematics()
{
////////////////////////////////////
- const Kinematics::Limits eta_limits = cuts_.central_cuts[Cuts::eta_single];
- y_min_ = eta_limits.min();
- y_max_ = eta_limits.max();
+ const Kinematics::Limits rap_limits = cuts_.central_cuts[Cuts::rapidity_single];
+ y_min_ = rap_limits.min();
+ y_max_ = rap_limits.max();
///////////// FIXME ////////////////
// Outgoing leptons
y1_ = y_min_+( y_max_-y_min_ )*xkt( 0 );
y2_ = y_min_+( y_max_-y_min_ )*xkt( 1 );
DebuggingInsideLoop( Form( "leptons rapidities (%.2f < y < %.2f): %f / %f", y_min_, y_max_, y1_, y2_ ) );
Kinematics::Limits ptdiff_limits = cuts_.central_cuts[Cuts::pt_diff];
if ( !ptdiff_limits.hasMax() ) ptdiff_limits.max() = 500.; //FIXME
pt_diff_ = ptdiff_limits.min()+( ptdiff_limits.range() )*xkt( 2 );
phi_pt_diff_ = 2.*M_PI*xkt( 3 );
DebuggingInsideLoop( Form( "leptons pt difference:\n\t"
" mag = %f (%.2f < Dpt < %.2f)\n\t"
" phi = %f",
pt_diff_, ptdiff_limits.min(), ptdiff_limits.max(), phi_pt_diff_ ) );
}
double
PPtoLL::computeJacobian()
{
double jac = GenericKTProcess::minimalJacobian();
jac *= ( y_max_-y_min_ ); // d(y1)
jac *= ( y_max_-y_min_ ); // d(y2)
jac *= cuts_.central_cuts[Cuts::pt_diff].range(); // d(Dpt)
jac *= 2.*M_PI; // d(phiDpt)
return jac;
}
double
PPtoLL::computeKTFactorisedMatrixElement()
{
const double mp = Particle::massFromPDGId( Particle::Proton ), mp2 = mp*mp;
const double ml = event_->getByRole( Particle::CentralSystem )[0].mass(), ml2 = ml*ml;
const unsigned int iterm11 = 1, // Long-long
iterm22 = 1, // Trans-trans
iterm12 = 1, // Long-trans
itermtt = 1; // Trans-trans(')
//=================================================================
// How matrix element is calculated
//=================================================================
const bool off_shell = true;
//=================================================================
// two terms in Wolfgang's formula for
// off-shell gamma gamma --> l^+ l^-
//=================================================================
const unsigned int imat1 = 2, imat2 = 0;
//=================================================================
// matrix element computation
//=================================================================
//const double stild = s_/2.*(1+sqrt(1.-(4*pow(mp2, 2))/s_*s_));
// Inner photons
const double q1tx = qt1_*cos( phi_qt1_ ), q1ty = qt1_*sin( phi_qt1_ ),
q2tx = qt2_*cos( phi_qt2_ ), q2ty = qt2_*sin( phi_qt2_ );
DebuggingInsideLoop( Form( "q1t(x/y) = %e / %e\n\t"
"q2t(x/y) = %e / %e", q1tx, q1ty, q2tx, q2ty ) );
// Two-photon system
const double ptsumx = q1tx+q2tx,
ptsumy = q1ty+q2ty,
ptsum = sqrt( ptsumx*ptsumx+ptsumy*ptsumy );
const double ptdiffx = pt_diff_*cos( phi_pt_diff_ ),
ptdiffy = pt_diff_*sin( phi_pt_diff_ );
// Outgoing leptons
const double pt1x = ( ptsumx+ptdiffx )*0.5, pt1y = ( ptsumy+ptdiffy )*0.5, pt1 = sqrt( pt1x*pt1x+pt1y*pt1y ),
pt2x = ( ptsumx-ptdiffx )*0.5, pt2y = ( ptsumy-ptdiffy )*0.5, pt2 = sqrt( pt2x*pt2x+pt2y*pt2y );
const Kinematics::Limits pt_limits = cuts_.central_cuts[Cuts::pt_single];
if ( pt_limits.hasMin() && ( pt1 < pt_limits.min() || pt2 < pt_limits.min() ) ) return 0.;
if ( pt_limits.hasMax() && ( pt1 > pt_limits.max() || pt2 > pt_limits.max() ) ) return 0.;
// transverse mass for the two leptons
const double amt1 = sqrt( pt1*pt1+ml2 ),
amt2 = sqrt( pt2*pt2+ml2 );
//=================================================================
// a window in transverse momentum difference
//=================================================================
const Kinematics::Limits ptdiff_limits = cuts_.central_cuts[Cuts::pt_diff];
if ( ptdiff_limits.hasMax() && fabs( pt1-pt2 ) > ptdiff_limits.max() ) return 0.;
//=================================================================
// a window in rapidity distance
//=================================================================
const double dely = fabs( y1_-y2_ );
const Kinematics::Limits dely_limits = cuts_.central_cuts[Cuts::dely];
if ( dely_limits.hasMin() && dely < dely_limits.min() ) return 0.;
if ( dely_limits.hasMax() && dely > dely_limits.max() ) return 0.;
//=================================================================
// auxiliary quantities
//=================================================================
const double alpha1 = amt1/sqs_*exp( y1_ ),
alpha2 = amt2/sqs_*exp( y2_ ),
beta1 = amt1/sqs_*exp( -y1_ ),
beta2 = amt2/sqs_*exp( -y2_ );
DebuggingInsideLoop( Form( "Sudakov parameters:\n\t"
" alpha1/2 = %f / %f\n\t"
" beta1/2 = %f / %f", alpha1, alpha2, beta1, beta2 ) );
const double q1t2 = q1tx*q1tx+q1ty*q1ty,
q2t2 = q2tx*q2tx+q2ty*q2ty;
//const double old_x2 = 0.; //FIXME figure out where this comes from
//const double delta_x1 = (MX_*MX_+q2t2)/((1.-old_x2)*s_);
//x1 = alpha1+alpha2+delta_x1;
const double x1 = alpha1+alpha2,
x2 = beta1 +beta2;
/*const double xi_x1 = log10(x1);
const double xi_x2 = log10(x2);*/
const double z1p = alpha1/x1, z1m = alpha2/x1,
z2p = beta1 /x2, z2m = beta2 /x2;
DebuggingInsideLoop( Form( "z(1/2)p = %f / %f\n\t"
"z(1/2)m = %f / %f", z1p, z2p, z1m, z2m ) );
if ( x1 > 1. || x2 > 1. ) return 0.; // sanity check
// FIXME FIXME FIXME
const double ak10 = event_->getOneByRole( Particle::IncomingBeam1 ).energy(),
ak1z = event_->getOneByRole( Particle::IncomingBeam1 ).momentum().pz(),
ak20 = event_->getOneByRole( Particle::IncomingBeam2 ).energy(),
ak2z = event_->getOneByRole( Particle::IncomingBeam2 ).momentum().pz();
DebuggingInsideLoop( Form( "incoming particles: p1: %f / %f\n\t"
" p2: %f / %f", ak1z, ak10, ak2z, ak20 ) );
//=================================================================
// additional conditions for energy-momentum conservation
//=================================================================
const double s1_eff = x1*s_-qt1_*qt1_, s2_eff = x2*s_-qt2_*qt2_;
const double invm = sqrt( amt1*amt1 + amt2*amt2 + 2.*amt1*amt2*cosh(y1_-y2_) - ptsum*ptsum );
DebuggingInsideLoop( Form( "s(1/2)_eff = %f / %f GeV^2\n\t"
"dilepton invariant mass = %f GeV", s1_eff, s2_eff, invm ) );
switch ( cuts_.mode ) {
case Kinematics::ElasticInelastic: if ( sqrt( s1_eff ) <= ( MY_+invm ) ) return 0.;
case Kinematics::InelasticElastic: if ( sqrt( s2_eff ) <= ( MX_+invm ) ) return 0.;
case Kinematics::InelasticInelastic: if ( sqrt( s1_eff ) <= ( MY_+invm ) ) return 0.;
if ( sqrt( s2_eff ) <= ( MX_+invm ) ) return 0.;
default: break;
}
//const double qcaptx = pcaptx, qcapty = pcapty;
//=================================================================
// four-momenta of the outgoing protons (or remnants)
//=================================================================
const double px_plus = ( 1.-x1 )*fabs( ak1z )*sqrt( 2. ),
px_minus = ( MX_*MX_ + q1tx*q1tx + q1ty*q1ty )*0.5/px_plus;
const double py_minus = ( 1.-x2 )*fabs( ak2z )*sqrt( 2. ), // warning! sign of pz??
py_plus = ( MY_*MY_ + q2tx*q2tx + q2ty*q2ty )*0.5/py_minus;
DebuggingInsideLoop( Form( "px_(+/-) = %f / %f\n\t"
"py_(+/-) = %f / %f", px_plus, px_minus, py_plus, py_minus ) );
PX_ = Particle::Momentum( -q1tx, -q1ty, ( px_plus-px_minus )*0.5*sqrt( 2. ), (px_plus+px_minus)*0.5*sqrt( 2. ) );
PY_ = Particle::Momentum( -q2tx, -q2ty, ( py_plus-py_minus )*0.5*sqrt( 2. ), (py_plus+py_minus)*0.5*sqrt( 2. ) );
DebuggingInsideLoop( Form( "First remnant: (E,p) = (%f, %f, %f, %f)\n\t"
"Second remnant: (E,p) = (%f, %f, %f, %f)",
PX_.px(), PX_.py(), PX_.pz(), PX_.energy(),
PY_.px(), PY_.py(), PY_.pz(), PY_.energy() ) );
assert( fabs( PX_.mass()-MX_ ) < 1.e-6 );
assert( fabs( PY_.mass()-MY_ ) < 1.e-6 );
//=================================================================
// four-momenta of the outgoing l^+ and l^-
//=================================================================
Particle::Momentum p1( pt1x, pt1y, alpha1*ak1z + beta1*ak2z, alpha1*ak10 + beta1*ak20 ),
p2( pt2x, pt2y, alpha2*ak1z + beta2*ak2z, alpha2*ak10 + beta2*ak20 );
DebuggingInsideLoop( Form( "unboosted first lepton: (E,p), m = (%f, %f, %f, %f), %f\n\t"
" second lepton: (E,p), m = (%f, %f, %f, %f), %f",
p1.px(), p1.py(), p1.pz(), p1.energy(), p1.mass(),
p2.px(), p2.py(), p2.pz(), p2.energy(), p2.mass() ) );
Pl1_ = Particle::Momentum( pt1x, pt1y, sqrt( pt1*pt1 + ml2 )*sinh( y1_ ), sqrt( pt1*pt1 + ml2 )*cosh( y1_ ) );
Pl2_ = Particle::Momentum( pt2x, pt2y, sqrt( pt2*pt2 + ml2 )*sinh( y2_ ), sqrt( pt2*pt2 + ml2 )*cosh( y2_ ) );
DebuggingInsideLoop( Form( "First lepton: (E,p), m = (%f, %f, %f, %f), %f\n\t"
"Second lepton: (E,p), m = (%f, %f, %f, %f), %f",
Pl1_.px(), Pl1_.py(), Pl1_.pz(), Pl1_.energy(), Pl1_.mass(),
Pl2_.px(), Pl2_.py(), Pl2_.pz(), Pl2_.energy(), Pl2_.mass() ) );
assert( fabs( Pl1_.mass()-event_->getByRole( Particle::CentralSystem )[0].mass() ) < 1.e-6 );
assert( fabs( Pl2_.mass()-event_->getByRole( Particle::CentralSystem )[1].mass() ) < 1.e-6 );
//=================================================================
// four-momenta squared of the virtual photons
//=================================================================
// FIXME FIXME FIXME /////////////////////
Particle::Momentum q1( q1tx, q1ty, 0., 0. ),
q2( q2tx, q2ty, 0., 0. );
//////////////////////////////////////////
DebuggingInsideLoop( Form( "First photon*: (E,p), m2 = (%f, %f, %f, %f), %e\n\t"
"Second photon*: (E,p), m2 = (%f, %f, %f, %f), %e",
q1.px(), q1.py(), q1.pz(), q1.energy(), q1.mass2(),
q2.px(), q2.py(), q2.pz(), q2.energy(), q2.mass2() ) );
//const double q12 = q1.mass2(), q22 = q2.mass2();
//=================================================================
// Mendelstam variables
//=================================================================
//const double shat = s_*x1*x2; // ishat = 1 (approximation)
//const double shat = ( q1+q2 ).mass2(); // ishat = 2 (exact formula)
const double that1 = ( q1-p1 ).mass2(), that2 = ( q2-p2 ).mass2(),
uhat1 = ( q1-p2 ).mass2(), uhat2 = ( q2-p1 ).mass2();
DebuggingInsideLoop( Form( "that(1/2) = %f / %f\n\t"
"uhat(1/2) = %f / %f",
that1, that2, uhat1, uhat2 ) );
//const double mll = sqrt( shat );
const double that = 0.5*( that1+that2 ), uhat = 0.5*( uhat1+uhat2 );
//=================================================================
// matrix elements
//=================================================================
double amat2 = 0.;
if ( !off_shell ) {
//=================================================================
// on-shell formula for M^2
//=================================================================
const double ml4 = ml2*ml2, ml8 = ml4*ml4;
const double term1 = 6. *ml8,
term2 = -3. *ml4*that*that,
term3 = -14.*ml4*that*uhat,
term4 = -3. *ml4*uhat*uhat,
term5 = ml2*that*that*that,
term6 = 7.* ml2*that*that*uhat,
term7 = 7.* ml2*that*uhat*uhat,
term8 = ml2*uhat*uhat*uhat,
term9 = -that*that*that*uhat,
term10 = -that*uhat*uhat*uhat;
const double auxil_gamgam = -2.*( term1+term2+term3+term4+term5+term6+term7+term8+term9+term10 )/( pow( ( ml2-that )*( ml2-uhat ), 2) );
const double g_em_sq = 4.*M_PI*Constants::alphaEM;
amat2 = g_em_sq*g_em_sq*auxil_gamgam;
}
else if ( off_shell ) {
//=================================================================
// Wolfgang's formulae
//=================================================================
double aux2_1, aux2_2;
const double ak1_x = z1m*pt1x - z1p*pt2x, ak1_y = z1m*pt1y - z1p*pt2y,
ak2_x = z2m*pt1x - z2p*pt2x, ak2_y = z2m*pt1y - z2p*pt2y;
const double t1abs = ( q1t2 + x1*( MX_*MX_-mp2 )+x1*x1*mp2 )/( 1.-x1 ),
t2abs = ( q2t2 + x2*( MY_*MY_-mp2 )+x2*x2*mp2 )/( 1.-x2 );
const double eps12 = ml2 + z1p*z1m*t1abs,
eps22 = ml2 + z2p*z2m*t2abs;
const double Phi10 = 1./( pow( ak1_x + z1p*q2tx, 2 ) + pow( ak1_y + z1p*q2ty, 2 ) + eps12 )
-1./( pow( ak1_x - z1m*q2tx, 2 ) + pow( ak1_y - z1m*q2ty, 2 ) + eps12 ),
Phi11_x = ( ak1_x + z1p*q2tx )/( pow( ak1_x + z1p*q2tx, 2 ) + pow( ak1_y + z1p*q2ty, 2 ) + eps12 )
-( ak1_x - z1m*q2tx )/( pow( ak1_x - z1m*q2tx, 2 ) + pow( ak1_y - z1m*q2ty, 2 ) + eps12 ),
Phi11_y = ( ak1_y + z1p*q2ty )/( pow( ak1_x + z1p*q2tx, 2 ) + pow( ak1_y + z1p*q2ty, 2 ) + eps12 )
-( ak1_y - z1m*q2ty )/( pow( ak1_x - z1m*q2tx, 2 ) + pow( ak1_y - z1m*q2ty, 2 ) + eps12 ),
Phi102 = Phi10*Phi10;
//const double Phi112 = Phi11_x*Phi11_x + Phi11_y*Phu11_y;
const double Phi20 = 1./( pow( ak2_x + z2p*q1tx, 2 )+pow( ak2_y + z2p*q1ty, 2 ) + eps22 )
-1./( pow( ak2_x - z2m*q1tx, 2 )+pow( ak2_y - z2m*q1ty, 2 ) + eps22 ),
Phi21_x = ( ak2_x + z2p*q1tx )/( pow( ak2_x + z2p*q1tx, 2 ) + pow( ak2_y + z2p*q1ty, 2 ) + eps22 )
-( ak2_x - z2m*q1tx )/( pow( ak2_x - z2m*q1tx ,2 ) + pow( ak2_y - z2m*q1ty, 2 ) + eps22 ),
Phi21_y = ( ak2_y + z2p*q1ty )/( pow( ak2_x + z2p*q1tx, 2 ) + pow( ak2_y + z2p*q1ty, 2 ) + eps22 )
-( ak2_y - z2m*q1ty )/( pow( ak2_x - z2m*q1tx, 2 ) + pow( ak2_y - z2m*q1ty, 2 ) + eps22 ),
Phi202 = Phi20*Phi20;
//const double Phi212 = Phi21_x*Phi21_x + Phi21_y*Phi21_y;
/*aux2_1 = iterm11*( ml2 + 4.*z1p*z1m*t1abs )*Phi102
+iterm22*( z1p*z1p + z1m*z1m )*Phi112
-iterm12*4.*z1p*z1m*( z1p-z1m )*Phi10*( q1tx*Phi11_x + q1ty*Phi11_y );
aux2_2 = iterm11*( ml2+4.*z2p*z2m*t2abs)*Phi202
+iterm22*( z2p*z2p + z2m*z2m )*Phi212
-iterm12*4.*z2p*z2m*( z2p-z2m )*Phi20*( q2tx*Phi21_x + q2ty*Phi21_y );*/
const double Phi11_dot_e = ( Phi11_x*q1tx + Phi11_y*q1ty )/qt1_, Phi11_cross_e = ( Phi11_x*q1ty-Phi11_y*q1tx )/qt1_;
const double Phi21_dot_e = ( Phi21_x*q2tx + Phi21_y*q2ty )/qt2_, Phi21_cross_e = ( Phi21_x*q2ty-Phi21_y*q2tx )/qt2_;
DebuggingInsideLoop( Form( "Phi1: E, px, py = %e, %e, %e\n\t"
"Phi2: E, px, py = %e, %e, %e\n\t"
"(dot): %e / %e\n\t"
"(cross): %e / %e",
Phi10, Phi11_x, Phi11_y, Phi20, Phi21_x, Phi21_y,
Phi11_dot_e, Phi21_dot_e, Phi11_cross_e, Phi21_cross_e ) );
aux2_1 = iterm11 * ( ml2 + 4.*z1p*z1p*z1m*z1m*t1abs ) * Phi102
+iterm22 * ( ( z1p*z1p + z1m*z1m )*( Phi11_dot_e*Phi11_dot_e + Phi11_cross_e*Phi11_cross_e ) )
+itermtt * ( Phi11_cross_e*Phi11_cross_e - Phi11_dot_e*Phi11_dot_e )
-iterm12 * 4.*z1p*z1m*( z1p-z1m ) * Phi10 * ( q1tx*Phi11_x + q1ty*Phi11_y );
aux2_2 = iterm11 * ( ml2 + 4.*z2p*z2p*z2p*z2m*t2abs ) * Phi202
+iterm22 * ( ( z2p*z2p + z2m*z2m )*( Phi21_dot_e*Phi21_dot_e + Phi21_cross_e*Phi21_cross_e ) )
+itermtt * ( Phi21_cross_e*Phi21_cross_e - Phi21_dot_e*Phi21_dot_e )
-iterm12 * 4.*z2p*z2m*( z2p-z2m ) * Phi20 * ( q2tx*Phi21_x + q2ty*Phi21_y );
//=================================================================
// convention of matrix element as in our kt-factorization
// for heavy flavours
//=================================================================
const double norm = 16.*M_PI*M_PI*Constants::alphaEM*Constants::alphaEM;
double amat2_1 = norm * pow( x1*x2*s_, 2 ) * aux2_1*2.*z1p*z1m*t1abs/( q1t2*q2t2 ) * t2abs/q2t2;
double amat2_2 = norm * pow( x1*x2*s_, 2 ) * aux2_2*2.*z2p*z2m*t2abs/( q1t2*q2t2 );
//=================================================================
// symmetrization
//=================================================================
amat2 = 0.5*( imat1*amat2_1 + imat2*amat2_2 );
DebuggingInsideLoop( Form( "aux2(1/2) = %e / %e\n\t"
"amat2(1/2), amat2 = %e / %e / %e", aux2_1, aux2_2, amat2_1, amat2_2, amat2 ) );
/*const double xx1 = alpha1+alpha2, xx2 = beta1+beta2;
const double sudakov_2 = ( MX_*MX_ - mp2+q2t2+xx2*mp2 )/( ( 1.-xx2 )*s_ );
const double sudakov_1 = ( q1t2 + xx1*mp2 )/( ( 1.-xx1 )*s_ );
const double ratio1 = sudakov_1 / xx1,
ratio2 = sudakov_2 / xx2;*/
//if ( ratio1 > 0.01 ) return 0.;
}
//============================================
// unintegrated photon distributions
//============================================
GenericKTProcess::computeIncomingFluxes( x1, q1t2, x2, q2t2 );
//=================================================================
// factor 2.*pi from integration over phi_sum
// factor 1/4 from jacobian of transformations
// factors 1/pi and 1/pi due to integration over
// d^2 kappa_1 d^2 kappa_2 instead d kappa_1^2 d kappa_2^2
//=================================================================
const double aintegral = amat2 / ( 16.*M_PI*M_PI*x1*x1*x2*x2*s_*s_ )
* flux1_/M_PI * flux2_/M_PI
* Constants::GeV2toBarn * 0.25;
//=================================================================
return aintegral*qt1_*qt2_*pt_diff_;
//=================================================================
}
void
PPtoLL::fillCentralParticlesKinematics()
{
// randomise the charge of the outgoing leptons
int sign = ( drand()>.5 ) ? +1 : -1;
//=================================================================
// first outgoing lepton
//=================================================================
Particle& ol1 = event_->getByRole( Particle::CentralSystem )[0];
ol1.setPdgId( ol1.pdgId(), sign );
ol1.setStatus( Particle::FinalState );
ol1.setMomentum( Pl1_ );
//=================================================================
// second outgoing lepton
//=================================================================
Particle& ol2 = event_->getByRole( Particle::CentralSystem )[1];
ol2.setPdgId( ol2.pdgId(), -sign );
ol2.setStatus( Particle::FinalState );
ol2.setMomentum( Pl2_ );
}
diff --git a/test/test_physics_processes.cpp b/test/test_physics_processes.cpp
index a7cbad8..72fed77 100644
--- a/test/test_physics_processes.cpp
+++ b/test/test_physics_processes.cpp
@@ -1,123 +1,124 @@
#include "CepGen/Generator.h"
#include "CepGen/Processes/GamGamLL.h"
#include "CepGen/Processes/PPtoLL.h"
#include <assert.h>
int
main( int argc, char* argv[] )
{
typedef std::map<std::string,std::pair<double,double> > KinematicsMap;
typedef std::map<float, KinematicsMap> ValuesAtCutMap;
// values defined at pt(single lepton)>15 GeV, |eta(single lepton)|<2.5, mX<1000 GeV
// process -> { pt cut -> { kinematics -> ( sigma, delta(sigma) ) } }
std::map<std::string,ValuesAtCutMap> values_map = {
//--- LPAIR values at sqrt(s) = 13 TeV
{ "1_lpair", {
{ 3.0, { // pt cut
{ "1_elastic", { 2.0871703e1, 3.542e-2 } },
{ "2_singlediss", { 1.5042536e1, 3.256e-2 } },
{ "3_doublediss", { 1.38835e1, 4.03018e-2 } }
} },
{ 15.0, { // pt cut
{ "1_elastic", { 4.1994803e-1, 8.328e-4 } },
{ "2_singlediss", { 4.8504819e-1, 1.171e-3 } },
{ "3_doublediss", { 6.35650e-1, 1.93968e-3 } }
} },
} },
//--- PPTOLL values
{ "2_pptoll", {
{ 3.0, { // pt cut
{ "1_elastic", { 2.0936541e1, 1.4096e-2 } },
{ "2_singlediss", { 1.4844881e1, 2.0723e-2 } }, // SU, qt<50
{ "3_doublediss", { 1.3580637e1, 2.2497e-2 } }, // SU, qt<50
} },
{ 15.0, { // pt cut
{ "1_elastic", { 4.2515888e-1, 3.0351e-4 } },
{ "2_singlediss", { 4.4903253e-1, 5.8970e-4 } }, // SU, qt<50
{ "3_doublediss", { 5.1923819e-1, 9.6549e-4 } }, // SU, qt<50
/*{ "2_singlediss", { 4.6710287e-1, 6.4842e-4 } }, // SU, qt<500
{ "3_doublediss", { 5.6316353e-1, 1.1829e-3 } }, // SU, qt<500*/
} },
} },
};
const double num_sigma = 3.0;
if ( argc == 1 || strcmp( argv[1], "debug" ) != 0 ) {
CepGen::Logger::get().level = CepGen::Logger::Nothing;
}
Timer tmr;
CepGen::Generator mg;
mg.parameters->kinematics.setSqrtS( 13.e3 );
mg.parameters->kinematics.central_cuts[CepGen::Cuts::eta_single].in( -2.5, 2.5 );
mg.parameters->kinematics.remnant_cuts[CepGen::Cuts::mass].max() = 1000.;
//mg.parameters->vegas.ncvg = 50000;
mg.parameters->vegas.itvg = 5;
Information( Form( "Initial configuration time: %.3f ms", tmr.elapsed()*1.e3 ) );
tmr.reset();
unsigned short num_tests = 0, num_tests_passed = 0;
std::vector<std::string> failed_tests;
for ( const auto& values_vs_generator : values_map ) { // loop over all generators
if ( values_vs_generator.first == "1_lpair" ) mg.parameters->setProcess( new CepGen::Process::GamGamLL );
else if ( values_vs_generator.first == "2_pptoll" ) {
mg.parameters->setProcess( new CepGen::Process::PPtoLL );
mg.parameters->kinematics.structure_functions = CepGen::StructureFunctionsType::SzczurekUleshchenko; //FIXME move to a dedicated class
mg.parameters->kinematics.initial_cuts[CepGen::Cuts::qt].max() = 50.0;
+ mg.parameters->kinematics.central_cuts[CepGen::Cuts::rapidity_single].in( -2.5, 2.5 );
}
else { InError( Form( "Unrecognized generator mode: %s", values_vs_generator.first.c_str() ) ); break; }
for ( const auto& values_vs_cut : values_vs_generator.second ) { // loop over the single lepton pT cut
mg.parameters->kinematics.central_cuts[CepGen::Cuts::pt_single].min() = values_vs_cut.first;
for ( const auto& values_vs_kin : values_vs_cut.second ) { // loop over all possible kinematics
if ( values_vs_kin.first == "1_elastic" ) mg.parameters->kinematics.mode = CepGen::Kinematics::ElasticElastic;
else if ( values_vs_kin.first == "2_singlediss" ) mg.parameters->kinematics.mode = CepGen::Kinematics::InelasticElastic;
else if ( values_vs_kin.first == "3_doublediss" ) mg.parameters->kinematics.mode = CepGen::Kinematics::InelasticInelastic;
else { InError( Form( "Unrecognized kinematics mode: %s", values_vs_kin.first.c_str() ) ); break; }
Information( Form( "Process: %s/%s\n\tConfiguration time: %.3f ms", values_vs_generator.first.c_str(), values_vs_kin.first.c_str(), tmr.elapsed()*1.e3 ) );
tmr.reset();
mg.clearRun();
const double xsec_ref = values_vs_kin.second.first, err_xsec_ref = values_vs_kin.second.second;
double xsec_cepgen, err_xsec_cepgen;
mg.computeXsection( xsec_cepgen, err_xsec_cepgen );
const double sigma = ( fabs( xsec_ref-xsec_cepgen ) ) / sqrt( err_xsec_cepgen*err_xsec_cepgen + err_xsec_ref*err_xsec_ref );
Information( Form( "Computed cross section:\n\tRef. = %.3e +/- %.3e\n\tCepGen = %.3e +/- %.3e\n\tPull: %.6f", xsec_ref, err_xsec_ref, xsec_cepgen, err_xsec_cepgen, sigma ) );
Information( Form( "Computation time: %.3f ms", tmr.elapsed()*1.e3 ) );
tmr.reset();
if ( fabs( sigma )<num_sigma ) num_tests_passed++;
else {
failed_tests.emplace_back( values_vs_generator.first+":"+
Form( "pt-gt-%.1f", values_vs_cut.first )+":"+
values_vs_kin.first+":"
"ref="+Form( "%.3e", xsec_ref )+":"
"got="+Form( "%.3e", xsec_cepgen )+":"
"pull="+Form( "%.3f", sigma ) );
}
num_tests++;
std::cout << "Test " << num_tests_passed << "/" << num_tests << " passed!" << std::endl;
}
}
}
if ( failed_tests.size() != 0 ) {
std::ostringstream os;
for ( const auto& fail : failed_tests ) os << " (*) " << fail << std::endl;
throw CepGen::Exception( __PRETTY_FUNCTION__, Form( "Some tests failed!\n%s", os.str().c_str() ), CepGen::FatalError );
}
Information( "ALL TESTS PASSED!" );
return 0;
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, May 3, 5:53 AM (17 h, 11 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4982854
Default Alt Text
(86 KB)

Event Timeline