Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F7878486
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
48 KB
Subscribers
None
View Options
diff --git a/FixedOrderGen/include/EventGenerator.hh b/FixedOrderGen/include/EventGenerator.hh
index 4399336..cc63bfb 100644
--- a/FixedOrderGen/include/EventGenerator.hh
+++ b/FixedOrderGen/include/EventGenerator.hh
@@ -1,47 +1,50 @@
#pragma once
#include "RHEJ/PDF.hh"
#include "RHEJ/MatrixElement.hh"
#include "JetParameters.hh"
#include "Process.hh"
#include "Beam.hh"
#include "Status.hh"
+#include "HiggsProperties.hh"
namespace RHEJ{
class Event;
class HiggsCouplingSettings;
class ScaleGenerator;
}
namespace HEJFOG{
class EventGenerator{
public:
EventGenerator(
Process process,
Beam beam,
RHEJ::ScaleGenerator & scale_gen,
JetParameters jets,
int pdf_id,
double uno_chance,
+ HiggsProperties higgs_properties,
RHEJ::HiggsCouplingSettings Higgs_coupling
);
RHEJ::Event gen_event();
Status status() const {
return status_;
}
private:
RHEJ::PDF pdf_;
RHEJ::MatrixElement ME_;
RHEJ::ScaleGenerator & scale_gen_;
Process process_;
JetParameters jets_;
Beam beam_;
Status status_;
double uno_chance_;
+ HiggsProperties higgs_properties_;
};
}
diff --git a/FixedOrderGen/include/Higgs_properties.hh b/FixedOrderGen/include/Higgs_properties.hh
new file mode 100644
index 0000000..07d580b
--- /dev/null
+++ b/FixedOrderGen/include/Higgs_properties.hh
@@ -0,0 +1,11 @@
+#pragma once
+
+#include <vector>
+
+namespace HEJFOG{
+ struct HiggsProperties{
+ double mass;
+ double width;
+ std::vector<Decay> decays;
+ }
+}
diff --git a/FixedOrderGen/include/PhaseSpacePoint.hh b/FixedOrderGen/include/PhaseSpacePoint.hh
index 5ce0038..7fa3d1e 100644
--- a/FixedOrderGen/include/PhaseSpacePoint.hh
+++ b/FixedOrderGen/include/PhaseSpacePoint.hh
@@ -1,131 +1,145 @@
/** \file PhaseSpacePoint.hh
* \brief Contains the PhaseSpacePoint Class
*/
#pragma once
#include <vector>
#include <CLHEP/Random/Randomize.h>
#include <CLHEP/Random/RanluxEngine.h>
#include "fastjet/JetDefinition.hh"
#include "RHEJ/utility.hh"
#include "RHEJ/Event.hh"
#include "RHEJ/PDG_codes.hh"
#include "RHEJ/PDF.hh"
#include "Status.hh"
+#include "HiggsProperties.hh"
namespace HEJFOG{
class Process;
using RHEJ::Sparticle;
//! A point in resummation phase space
class PhaseSpacePoint{
public:
//! Default PhaseSpacePoint Constructor
PhaseSpacePoint() = default;
//! PhaseSpacePoint Constructor
/**
* @param proc The process to generate
* @param jet_def The Jet Definition Used
* @param jetptmin Minimum Jet Transverse Momentum
* @param rapmax Maximum parton rapidity
* @param pdf The pdf set (used for sampling)
* @param uno_chance Chance to turn a potentially unordered
* emission into an actual one
*
* Initially, only FKL phase space points are generated. If the most
* extremal emission in any direction is a quark or anti-quark and the
* next emission is a gluon, uno_chance is the chance to turn this into
* an unordered emission event by swapping the two flavours. At most one
* unordered emission will be generated in this way.
*/
PhaseSpacePoint(
Process const & proc,
fastjet::JetDefinition jet_def,double jetptmin, double rapmax,
RHEJ::PDF & pdf, double E_beam,
- double uno_chance
+ double uno_chance,
+ HiggsProperties const & higgs_properties
);
//! Get Weight Function
/**
* @returns Weight of Event
*/
double weight() const{
return weight_;
}
Status status() const{
return status_;
}
//! Get Incoming Function
/**
* @returns Incoming Particles
*/
std::array<Sparticle, 2> const & incoming() const{
return incoming_;
}
//! Get Outgoing Function
/**
* @returns Outgoing Particles
*/
std::vector<Sparticle> const & outgoing() const{
return outgoing_;
}
+ std::unordered_map<int, std::vector<Sparticle>> const & decays() const{
+ return decays_;
+ }
+
static void reset_ranlux(std::string const & init_file);
static void reset_ranlux(char const * init_file);
private:
std::vector<fastjet::PseudoJet> gen_LO_partons(
int count, double ptmin, double ptmax, double rapmax
);
Sparticle gen_boson(
- RHEJ::ParticleID bosonid
+ RHEJ::ParticleID bosonid, double mass, double width
);
bool jets_ok(
std::vector<fastjet::PseudoJet> const & Born_jets,
std::vector<fastjet::PseudoJet> const & partons
) const;
void reconstruct_incoming(
std::array<RHEJ::ParticleID, 2> const & ids,
RHEJ::PDF & pdf, double E_beam
);
RHEJ::ParticleID generate_incoming_id(double x, double uf, RHEJ::PDF & pdf);
bool momentum_conserved(double ep) const;
RHEJ::Sparticle const & most_backward_FKL(
std::vector<RHEJ::Sparticle> const & partons
) const;
RHEJ::Sparticle const & most_forward_FKL(
std::vector<RHEJ::Sparticle> const & partons
) const;
RHEJ::Sparticle & most_backward_FKL(std::vector<RHEJ::Sparticle> & partons) const;
RHEJ::Sparticle & most_forward_FKL(std::vector<RHEJ::Sparticle> & partons) const;
bool extremal_FKL_ok(
std::vector<fastjet::PseudoJet> const & partons
) const;
double random_normal(double stddev);
void maybe_turn_to_uno(double chance);
+ std::vector<Sparticle> decay_boson(
+ RHEJ::Sparticle const & parent,
+ std::vector<Decay> const & decays
+ );
+ Decay select_decay_channel(std::vector<Decay> const & decays);
+
double weight_;
Status status_;
double jetptmin_;
fastjet::JetDefinition jet_def_;
std::array<Sparticle, 2> incoming_;
std::vector<Sparticle> outgoing_;
+ //! Particle decays in the format {outgoing index, decay products}
+ std::unordered_map<int, std::vector<Sparticle>> decays_;
static CLHEP::Ranlux64Engine ran_;
};
RHEJ::UnclusteredEvent to_UnclusteredEvent(PhaseSpacePoint const & psp);
}
diff --git a/FixedOrderGen/include/config.hh b/FixedOrderGen/include/config.hh
index c1c57c1..9775f0a 100644
--- a/FixedOrderGen/include/config.hh
+++ b/FixedOrderGen/include/config.hh
@@ -1,36 +1,36 @@
#pragma once
#include "yaml-cpp/yaml.h"
#include "fastjet/JetDefinition.hh"
#include "RHEJ/HiggsCouplingSettings.hh"
#include "RHEJ/optional.hh"
#include "RHEJ/config.hh"
#include "RHEJ/output_formats.hh"
#include "RHEJ/exceptions.hh"
#include "Process.hh"
#include "JetParameters.hh"
#include "Beam.hh"
#include "HiggsProperties.hh"
namespace HEJFOG{
struct Config{
Process process;
int trials;
JetParameters jets;
Beam beam;
int pdf_id;
double unordered_fraction;
- HiggsProperties higgs_properties;
+ HiggsProperties Higgs_properties;
YAML::Node analysis_parameters;
RHEJ::ScaleGenerator scale_gen;
std::vector<RHEJ::OutputFile> output;
RHEJ::optional<std::string> RanLux_init;
RHEJ::HiggsCouplingSettings Higgs_coupling;
};
Config load_config(std::string const & config_file);
}
diff --git a/FixedOrderGen/src/EventGenerator.cc b/FixedOrderGen/src/EventGenerator.cc
index b72d377..6501f8d 100644
--- a/FixedOrderGen/src/EventGenerator.cc
+++ b/FixedOrderGen/src/EventGenerator.cc
@@ -1,73 +1,76 @@
#include "EventGenerator.hh"
#include "Process.hh"
#include "Beam.hh"
#include "JetParameters.hh"
#include "PhaseSpacePoint.hh"
#include "RHEJ/Event.hh"
#include "RHEJ/config.hh"
namespace HEJFOG{
EventGenerator::EventGenerator(
Process process,
Beam beam,
RHEJ::ScaleGenerator & scale_gen,
JetParameters jets,
int pdf_id,
double uno_chance,
+ HiggsProperties higgs_properties,
RHEJ::HiggsCouplingSettings Higgs_coupling
):
pdf_{pdf_id, beam.particles[0], beam.particles[1]},
ME_{
[this](double mu){ return pdf_.Halphas(mu); },
jets.def, jets.min_pt,
false,
std::move(Higgs_coupling)
},
scale_gen_{scale_gen},
process_{std::move(process)},
jets_{std::move(jets)},
beam_{std::move(beam)},
- uno_chance_{uno_chance}
+ uno_chance_{uno_chance},
+ higgs_properties_{std::move(higgs_properties)}
{
}
RHEJ::Event EventGenerator::gen_event(){
HEJFOG::PhaseSpacePoint psp{
process_,
jets_.def, jets_.min_pt, jets_.max_y,
- pdf_, beam_.energy,
- uno_chance_
+ pdf_, beam_.energy,
+ uno_chance_,
+ higgs_properties_
};
status_ = psp.status();
if(status_ != good) return {};
RHEJ::Event ev = scale_gen_(
RHEJ::Event{
to_UnclusteredEvent(std::move(psp)),
jets_.def, jets_.min_pt
}
);
const double shat = RHEJ::shat(ev);
const double xa = (ev.incoming()[0].E()-ev.incoming()[0].pz())/(2.*beam_.energy);
const double xb = (ev.incoming()[1].E()+ev.incoming()[1].pz())/(2.*beam_.energy);
// evaluate matrix element on this point
ev.central().weight *= ME_.tree(
ev.central().mur, ev.incoming(), ev.outgoing(), false
)/(shat*shat);
ev.central().weight *= pdf_.pdfpt(0,xa,ev.central().muf, ev.incoming()[0].type);
ev.central().weight *= pdf_.pdfpt(0,xb,ev.central().muf, ev.incoming()[1].type);
for(auto & var: ev.variations()){
var.weight *= ME_.tree(
var.mur, ev.incoming(), ev.outgoing(), false
)/(shat*shat);
var.weight *= pdf_.pdfpt(0,xa,var.muf, ev.incoming()[0].type);
var.weight *= pdf_.pdfpt(0,xb,var.muf, ev.incoming()[1].type);
}
return ev;
}
}
diff --git a/FixedOrderGen/src/PhaseSpacePoint.cc b/FixedOrderGen/src/PhaseSpacePoint.cc
index 716d498..6d74926 100644
--- a/FixedOrderGen/src/PhaseSpacePoint.cc
+++ b/FixedOrderGen/src/PhaseSpacePoint.cc
@@ -1,395 +1,447 @@
#include "PhaseSpacePoint.hh"
#include <random>
#include "RHEJ/kinematics.hh"
#include "RHEJ/utility.hh"
#include "RHEJ/debug.hh"
+#include "RHEJ/exceptions.hh"
#include "Process.hh"
using namespace RHEJ;
namespace HEJFOG{
static_assert(
std::numeric_limits<double>::has_quiet_NaN,
"no quiet NaN for double"
);
constexpr double NaN = std::numeric_limits<double>::quiet_NaN();
RHEJ::UnclusteredEvent to_UnclusteredEvent(PhaseSpacePoint const & psp){
RHEJ::UnclusteredEvent result;
result.incoming = psp.incoming();
std::sort(
begin(result.incoming), end(result.incoming),
[](Sparticle o1, Sparticle o2){return o1.p.pz()<o2.p.pz();}
);
assert(result.incoming.size() == 2);
result.outgoing = psp.outgoing();
assert(
std::is_sorted(
begin(result.outgoing), end(result.outgoing),
RHEJ::rapidity_less{}
)
);
assert(result.outgoing.size() >= 2);
+ result.decays = psp.decays();
result.central.mur = NaN;
result.central.muf = NaN;
result.central.weight = psp.weight();
return result;
}
namespace{
//generate Ranlux64Engine with fixed, predefined state
/*
* some (all?) of the Ranlux64Engine constructors leave fields
* uninitialised, invoking undefined behaviour. This can be
* circumvented by restoring the state from a file
*/
CLHEP::Ranlux64Engine gen_Ranlux64Engine(){
static const std::string state =
"9876\n"
"0.91280703978419097666\n"
"0.41606065829518357191\n"
"0.99156342622341142601\n"
"0.030922955274050423213\n"
"0.16206278421638486975\n"
"0.76151768001958330956\n"
"0.43765760066092695979\n"
"0.42904698253748563275\n"
"0.11476317525663759511\n"
"0.026620053590963976831\n"
"0.65953715764414511114\n"
"0.30136722624439826745\n"
"3.5527136788005009294e-15 4\n"
"1 202\n";
const std::string file = std::tmpnam(nullptr);
{
std::ofstream out{file};
out << state;
}
CLHEP::Ranlux64Engine result;
result.restoreStatus(file.c_str());
return result;
}
}
CLHEP::Ranlux64Engine PhaseSpacePoint::ran_{gen_Ranlux64Engine()};
void PhaseSpacePoint::reset_ranlux(std::string const & init_file){
reset_ranlux(init_file.c_str());
}
void PhaseSpacePoint::reset_ranlux(char const * init_file){
ran_.restoreStatus(init_file);
}
namespace{
bool can_swap_to_uno(
RHEJ::Sparticle const & p1, RHEJ::Sparticle const & p2
){
return is_parton(p1)
&& p1.type != RHEJ::pid::gluon
&& p2.type == RHEJ::pid::gluon;
}
}
void PhaseSpacePoint::maybe_turn_to_uno(double chance){
assert(outgoing_.size() >= 2);
const size_t nout = outgoing_.size();
const bool can_be_uno_backward = can_swap_to_uno(
outgoing_[0], outgoing_[1]
);
const bool can_be_uno_forward = can_swap_to_uno(
outgoing_[nout-1], outgoing_[nout-2]
);
if(!can_be_uno_backward && !can_be_uno_forward) return;
if(ran_.flat() < chance){
weight_ /= chance;
if(can_be_uno_backward && can_be_uno_forward){
if(ran_.flat() < 0.5){
std::swap(outgoing_[0].type, outgoing_[1].type);
}
else{
std::swap(outgoing_[nout-1].type, outgoing_[nout-2].type);
}
weight_ *= 2.;
}
else if(can_be_uno_backward){
std::swap(outgoing_[0].type, outgoing_[1].type);
}
else{
assert(can_be_uno_forward);
std::swap(outgoing_[nout-1].type, outgoing_[nout-2].type);
}
}
else weight_ /= 1 - chance;
}
PhaseSpacePoint::PhaseSpacePoint(
Process const & proc,
fastjet::JetDefinition jet_def,double jetptmin, double rapmax,
RHEJ::PDF & pdf, double E_beam,
- double uno_chance
+ double uno_chance,
+ HiggsProperties const & h
):
jetptmin_{jetptmin},
jet_def_{jet_def}
{
assert(proc.njets >= 2);
const int nout = proc.njets + (proc.boson?1:0);
status_ = good;
weight_ = 1;
weight_ /= std::tgamma(nout);
outgoing_.reserve(nout);
for(auto&& p_out: gen_LO_partons(proc.njets, jetptmin_, E_beam, rapmax)){
outgoing_.emplace_back(Sparticle{pid::gluon, std::move(p_out)});
}
if(status_ != good) return;
if(proc.boson && *proc.boson == pid::Higgs){
// The Higgs
- auto Hparticle=gen_boson(pid::higgs);
+ auto Hparticle=gen_boson(pid::higgs, h.mass, h.width);
auto pos=std::upper_bound(
begin(outgoing_),end(outgoing_),Hparticle,rapidity_less{}
);
- outgoing_.insert(pos,Hparticle);
+ outgoing_.insert(pos, Hparticle);
+ if(! h.decays.empty()){
+ const int boson_idx = std::distance(begin(outgoing_), pos);
+ decays_.emplace(boson_idx, decay_boson(outgoing_[boson_idx], h.decays));
+ }
}
reconstruct_incoming(proc.incoming, pdf, E_beam);
if(status_ != good) return;
// set outgoing states
most_backward_FKL(outgoing_).type = incoming_[0].type;
most_forward_FKL(outgoing_).type = incoming_[1].type;
if(proc.njets > 2) maybe_turn_to_uno(uno_chance);
}
std::vector<fastjet::PseudoJet> PhaseSpacePoint::gen_LO_partons(
int np , double ptmin, double ptmax, double rapmax
){
if (np<2) throw std::invalid_argument{"Not enough partons in gen_LO_partons"};
// heuristic parameters for pt sampling
const double ptpar = ptmin + np/5.;
const double arg_small_y = atan((ptmax - ptmin)/ptpar);
const double y_cut = 3.;
weight_ /= pow(16.*pow(M_PI,3),np-2);
std::vector<fastjet::PseudoJet> partons;
partons.reserve(np);
for(int i = 0; i < np; ++i){
const double y = -rapmax + 2*rapmax*ran_.flat();
weight_ *= 2*rapmax;
const double phi = 2*M_PI*ran_.flat();
weight_ *= 2.0*M_PI;
double pt;
const double r1 = ran_.flat();
if(y < y_cut){
pt = ptmin + ptpar*tan(r1*arg_small_y);
const double temp = cos(r1*arg_small_y);
weight_ *= pt*ptpar*arg_small_y/(temp*temp);
}
else{
const double ptpar2 = ptpar/(1 + 5*(y-y_cut));
const double temp = 1. - std::exp((ptmin-ptmax)/ptpar2);
pt = ptmin - ptpar2*std::log(1-r1*temp);
weight_ *= pt*ptpar2*temp/(1-r1*temp);
}
partons.emplace_back(fastjet::PtYPhiM(pt, y, phi));
assert(ptmin-1e-5 <= partons[i].pt() && partons[i].pt() <= ptmax+1e-5);
}
// Need to check that at LO, the number of jets = number of partons;
fastjet::ClusterSequence cs(partons, jet_def_);
auto cluster_jets=cs.inclusive_jets(jetptmin_);
if (cluster_jets.size()!=unsigned(np)){
weight_=0.0;
status_ = not_enough_jets;
return {};
}
std::sort(begin(partons), end(partons), rapidity_less{});
return partons;
}
Sparticle PhaseSpacePoint::gen_boson(
- RHEJ::ParticleID bosonid
+ RHEJ::ParticleID bosonid, double mass, double width
){
std::array<double,2> ptrans{0.,0.};
for (auto const & parton:outgoing_) {
ptrans[0]-=parton.px();
ptrans[1]-=parton.py();
}
// The Higgs:
// Generate a y Gaussian distributed around 0
const double y = random_normal(1.6);
- // r1=ran.flat();
- // double sH=flags.mH*(flags.mH + flags.GammaH*tan((M_PI*r1)/2. + (-1. + r1)*atan(flags.mH/flags.GammaH)));
- const double sH=125.*125.;
+ const double r1 = ran_.flat();
+ const double sH = mass*(
+ mass + width*tan(M_PI/2.*r1 + (r1-1.)*atan(mass/width))
+ );
const double mHperp=sqrt(ptrans[0]*ptrans[0]+ptrans[1]*ptrans[1]+sH);
const double pz=mHperp*sinh(y);
const double E=mHperp*cosh(y);
return Sparticle{bosonid,fastjet::PseudoJet{ptrans[0],ptrans[1],pz,E}};
}
Sparticle const & PhaseSpacePoint::most_backward_FKL(
std::vector<Sparticle> const & partons
) const{
if(!RHEJ::is_parton(partons[0])) return partons[1];
return partons[0];
}
Sparticle const & PhaseSpacePoint::most_forward_FKL(
std::vector<Sparticle> const & partons
) const{
const size_t last_idx = partons.size() - 1;
if(!RHEJ::is_parton(partons[last_idx])) return partons[last_idx-1];
return partons[last_idx];
}
Sparticle & PhaseSpacePoint::most_backward_FKL(
std::vector<Sparticle> & partons
) const{
if(!RHEJ::is_parton(partons[0])) return partons[1];
return partons[0];
}
Sparticle & PhaseSpacePoint::most_forward_FKL(
std::vector<Sparticle> & partons
) const{
const size_t last_idx = partons.size() - 1;
if(!RHEJ::is_parton(partons[last_idx])) return partons[last_idx-1];
return partons[last_idx];
}
void PhaseSpacePoint::reconstruct_incoming(
std::array<RHEJ::ParticleID, 2> const & ids,
RHEJ::PDF & pdf, double E_beam
){
std::tie(incoming_[0].p, incoming_[1].p) = incoming_momenta(outgoing_);
// calculate xa, xb
const double sqrts=2*E_beam;
const double xa=(incoming_[0].p.e()-incoming_[0].p.pz())/sqrts;
const double xb=(incoming_[1].p.e()+incoming_[1].p.pz())/sqrts;
// abort if phase space point is outside of collider energy reach
if (xa>1. || xb>1.){
weight_=0;
status_ = too_much_energy;
return;
}
// pick pdfs
/** TODO:
* ufa, ufb don't correspond to our final scale choice.
* The reversed HEJ scale generators currently expect a full event as input,
* so fixing this is not completely trivial
*/
if(ids[0] == pid::proton || ids[0] == pid::p_bar){
const double ufa=jetptmin_;
incoming_[0].type = generate_incoming_id(xa, ufa, pdf);
}
else {
incoming_[0].type = ids[0];
}
if(ids[1] == pid::proton || ids[1] == pid::p_bar){
const double ufb=jetptmin_;
incoming_[1].type = generate_incoming_id(xb, ufb, pdf);
}
else {
incoming_[1].type = ids[1];
}
assert(momentum_conserved(1e-7));
}
RHEJ::ParticleID PhaseSpacePoint::generate_incoming_id(
double x, double uf, RHEJ::PDF & pdf
){
const double pdfg=fabs(pdf.pdfpt(0,x,uf,pid::gluon));
const double pdfu=fabs(pdf.pdfpt(0,x,uf,pid::up));
const double pdfd=fabs(pdf.pdfpt(0,x,uf,pid::down));
const double pdfux=fabs(pdf.pdfpt(0,x,uf,pid::u_bar));
const double pdfdx=fabs(pdf.pdfpt(0,x,uf,pid::d_bar));
const double pdfc=fabs(pdf.pdfpt(0,x,uf,pid::charm));
const double pdfs=fabs(pdf.pdfpt(0,x,uf,pid::strange));
const double pdfsx=fabs(pdf.pdfpt(0,x,uf,pid::s_bar));
const double pdfb=fabs(pdf.pdfpt(0,x,uf,pid::b));
const double pdftot=pdfg+4.0/9.0*(pdfu + pdfd + pdfux + pdfdx +pdfs +pdfsx + 2.0*(pdfc +pdfb ));
const double r1=pdftot*ran_.flat();
double sum;
if (r1<(sum=pdfg)) {
weight_*=pdftot/pdfg;
return pid::gluon;
}
if (r1<(sum+=(4./9.)*pdfu)) {
weight_*=pdftot/((4./9.)*pdfu);
return pid::up;
}
else if (r1<(sum+=(4./9.)*pdfd)) {
weight_*=pdftot/((4./9.)*pdfd);
return pid::down;
}
else if (r1<(sum+=(4./9.)*pdfux)) {
weight_*=pdftot/((4./9.)*pdfux);
return pid::u_bar;
}
else if (r1<(sum+=(4./9.)*pdfdx)) {
weight_*=pdftot/((4./9.)*pdfdx);
return pid::d_bar;
}
else if (r1<(sum+=(4./9.)*pdfc)) {
weight_*=pdftot/((4./9.)*pdfc);
return pid::c;
}
else if (r1<(sum+=(4./9.)*pdfc)){
weight_*=pdftot/((4./9.)*pdfc);
return pid::c_bar;
}
else if (r1<(sum+=(4./9.)*pdfs)) {
weight_*=pdftot/((4./9.)*pdfs);
return pid::s;
}
else if (r1<(sum+=(4./9.)*pdfsx)) {
weight_*=pdftot/((4./9.)*pdfsx);
return pid::s_bar;
}
else if (r1<(sum+=(4./9.)*pdfb)) {
weight_*=pdftot/((4./9.)*pdfb);
return pid::b;
}
else if (r1<=(sum+=(4./9.)*pdfb)) {
weight_*=pdftot/((4./9.)*pdfb);
return pid::b_bar;
}
std::cout << "Error in choosing incoming parton: "<<x<<" "<<uf<<" "<<sum<<" "<<pdftot<<" "<<r1;
std::cout << " "<<pdfg+4./9.*(pdfu+pdfux+pdfd+pdfdx+pdfs+pdfsx+2.*(pdfc+pdfb))<<std::endl;
throw std::logic_error{"Failed to choose parton flavour"};
}
double PhaseSpacePoint::random_normal(double stddev){
const double r1 = ran_.flat();
const double r2 = ran_.flat();
const double lninvr1 = -log(r1);
const double result = stddev*sqrt(2.*lninvr1)*cos(2.*M_PI*r2);
weight_ *= exp(result*result/(2*stddev*stddev))*sqrt(2.*M_PI)*stddev;
return result;
}
bool PhaseSpacePoint::momentum_conserved(double ep) const{
fastjet::PseudoJet diff;
for(auto const & in: incoming()) diff += in.p;
for(auto const & out: outgoing()) diff -= out.p;
return nearby_ep(diff, fastjet::PseudoJet{}, ep);
}
+ Decay PhaseSpacePoint::select_decay_channel(std::vector<Decay> const & decays){
+ double br_total = 0.;
+ for(auto const & decay: decays) br_total += decay.branching_ratio;
+ // adjust weight
+ // this is given by (channel branching ratio)/(chance to pick channel)
+ // where (chance to pick channel) =
+ // (channel branching ratio)/(total branching ratio)
+ weight_ *= br_total;
+ const double r1 = br_total*ran_.flat();
+ double br_sum = 0.;
+ for(auto const & decay: decays){
+ br_sum += decay.branching_ratio;
+ if(r1 < br_sum) return decay;
+ }
+ throw std::logic_error{"unreachable"};
+ }
+
+ std::vector<Sparticle> PhaseSpacePoint::decay_boson(
+ RHEJ::Sparticle const & parent,
+ std::vector<Decay> const & decays
+ ){
+ const auto channel = select_decay_channel(decays);
+ if(channel.products.size() != 2){
+ throw RHEJ::not_implemented{
+ "only decays into two particles are implemented"
+ };
+ }
+ std::vector<Sparticle> decay_products(channel.products.size());
+ for(size_t i = 0; i < channel.products.size(); ++i){
+ decay_products[i].type = channel.products[i];
+ }
+ // choose polar and azimuth angle in parent rest frame
+ const double E = parent.m()/2;
+ const double theta = 2.*M_PI*ran_.flat();
+ const double cos_phi = 2.*ran_.flat()-1.;
+ const double sin_phi = sqrt(1. - cos_phi*cos_phi); // Know 0 < phi < pi
+ const double px = E*cos(theta)*sin_phi;
+ const double py = E*sin(theta)*sin_phi;
+ const double pz = E*cos_phi;
+ decay_products[0].p.reset(px, py, pz, E);
+ decay_products[1].p.reset(-px, -py, -pz, E);
+ for(auto & particle: decay_products) particle.p.boost(parent.p);
+ return decay_products;
+ }
}
diff --git a/FixedOrderGen/src/config.cc b/FixedOrderGen/src/config.cc
index 1454d63..cdd7d38 100644
--- a/FixedOrderGen/src/config.cc
+++ b/FixedOrderGen/src/config.cc
@@ -1,264 +1,264 @@
#include "config.hh"
#include <cctype>
#include "RHEJ/config.hh"
namespace HEJFOG{
using RHEJ::set_from_yaml;
using RHEJ::set_from_yaml_if_defined;
namespace{
//! Get YAML tree of supported options
/**
* The configuration file is checked against this tree of options
* in assert_all_options_known.
*/
YAML::Node const & get_supported_options(){
const static YAML::Node supported = [](){
YAML::Node supported;
static const auto opts = {
"process", "trials", "unordered fraction", "scales", "scale factors",
"max scale ratio", "pdf", "event output", "analysis", "RanLux init",
};
// add subnodes to "supported" - the assigned value is irrelevant
for(auto && opt: opts) supported[opt] = "";
for(auto && jet_opt: {"min pt", "algorithm", "R", "max rapidity"}){
supported["jets"][jet_opt] = "";
}
for(auto && Higgs_opt: {"mass", "width"}){
supported["Higgs properties"][Higgs_opt] = "";
}
supported["Higgs properties"]["decays"]["into"] = "";
supported["Higgs properties"]["decays"]["branching ratio"] = "";
for(auto && opt: {"mt", "use impact factors", "include bottom", "mb"}){
supported["Higgs coupling"][opt] = "";
}
for(auto && beam_opt: {"energy", "particles"}){
supported["beam"][beam_opt] = "";
}
return supported;
}();
return supported;
}
JetParameters get_jet_parameters(
YAML::Node const & node, std::string const & entry
){
const auto p = RHEJ::get_jet_parameters(node, entry);
JetParameters result;
result.def = p.def;
result.min_pt = p.min_pt;
set_from_yaml(result.max_y, node, entry, "max rapidity");
return result;
}
Beam get_Beam(
YAML::Node const & node, std::string const & entry
){
Beam beam;
std::vector<RHEJ::ParticleID> particles;
set_from_yaml(beam.energy, node, entry, "energy");
set_from_yaml_if_defined(particles, node, entry, "particles");
if(! particles.empty()){
for(RHEJ::ParticleID particle: particles){
if(particle != RHEJ::pid::p && particle != RHEJ::pid::p_bar){
throw std::invalid_argument{
"Unsupported value in option " + entry + ": particles:"
" only proton ('p') and antiproton ('p_bar') beams are supported"
};
}
}
if(particles.size() != 2){
throw std::invalid_argument{"Not exactly two beam particles"};
}
beam.particles.front() = particles.front();
beam.particles.back() = particles.back();
}
return beam;
}
std::vector<std::string> split(
std::string const & str, std::string const & delims
){
std::vector<std::string> result;
for(size_t begin, end = 0; end != str.npos;){
begin = str.find_first_not_of(delims, end);
if(begin == str.npos) break;
end = str.find_first_of(delims, begin + 1);
result.emplace_back(str.substr(begin, end - begin));
}
return result;
}
std::invalid_argument invalid_incoming(std::string const & what){
return std::invalid_argument{
"Incoming particle type " + what + " not supported,"
" incoming particles have to be 'p', 'p_bar' or partons"
};
}
std::invalid_argument invalid_outgoing(std::string const & what){
return std::invalid_argument{
"Outgoing particle type " + what + " not supported,"
" outgoing particles have to be 'j', 'photon', 'W+', 'W-', 'Z', 'H'"
};
}
Process get_process(
YAML::Node const & node, std::string const & entry
){
Process result;
std::string process_string;
set_from_yaml(process_string, node, entry);
assert(! process_string.empty());
const auto particles = split(process_string, " \n\t\v=>");
if(particles.size() < 3){
throw std::invalid_argument{
"Bad format in option process: '" + process_string
+ "', expected format is 'in1 in2 => out1 ...'"
};
}
result.incoming.front() = RHEJ::to_ParticleID(particles[0]);
result.incoming.back() = RHEJ::to_ParticleID(particles[1]);
for(size_t i = 0; i < result.incoming.size(); ++i){
const RHEJ::ParticleID in = result.incoming[i];
if(
in != RHEJ::pid::proton && in != RHEJ::pid::p_bar
&& !RHEJ::is_parton(in)
){
throw invalid_incoming(particles[i]);
}
}
result.njets = 0;
for(size_t i = result.incoming.size(); i < particles.size(); ++i){
assert(! particles[i].empty());
if(particles[i] == "j") ++result.njets;
else if(std::isdigit(particles[i].front())){
if(particles[i].back() != 'j'){
throw invalid_outgoing(particles[i]);
}
result.njets += std::stoi(particles[i]);
}
else{
const auto pid = RHEJ::to_ParticleID(particles[i]);
if(!RHEJ::is_AWZH_boson(pid)){
throw invalid_outgoing(particles[i]);
}
if(result.boson){
throw std::invalid_argument{
"More than one outgoing boson is not supported"
};
}
result.boson = pid;
}
}
if(result.njets < 2){
throw std::invalid_argument{
"Process has to include at least two jets ('j')"
};
}
return result;
}
Decay get_decay(YAML::Node const & node){
Decay decay;
set_from_yaml(decay.products, node, "into");
set_from_yaml(decay.branching_ratio, node, "branching ratio");
return decay;
}
std::vector<Decay> get_decays(YAML::Node const & node){
using YAML::NodeType;
if(!node) return {};
switch(node.Type()){
case NodeType::Null:
case NodeType::Undefined:
return {};
case NodeType::Scalar:
throw RHEJ::invalid_type{"value is not a list of decays"};
case NodeType::Map:
return {get_decay(node)};
case NodeType::Sequence:
std::vector<Decay> result;
for(auto && decay_str: node){
result.emplace_back();
set_from_yaml(result.back().products, decay_str, "into");
set_from_yaml(result.back().branching_ratio, decay_str, "branching ratio");
}
return result;
}
throw std::logic_error{"unreachable"};
}
HiggsProperties get_higgs_properties(
YAML::Node const & node, std::string const & entry
){
HiggsProperties result;
set_from_yaml(result.mass, node, entry, "mass");
set_from_yaml(result.width, node, entry, "width");
try{
result.decays = get_decays(node[entry]["decays"]);
}
catch(RHEJ::missing_option const & ex){
throw RHEJ::missing_option{entry + ": decays: " + ex.what()};
}
catch(RHEJ::invalid_type const & ex){
throw RHEJ::invalid_type{entry + ": decays: " + ex.what()};
}
return result;
}
Config to_Config(YAML::Node const & yaml){
try{
RHEJ::assert_all_options_known(yaml, get_supported_options());
}
catch(RHEJ::unknown_option const & ex){
throw RHEJ::unknown_option{std::string{"Unknown option '"} + ex.what() + "'"};
}
Config config;
config.process = get_process(yaml, "process");
set_from_yaml(config.trials, yaml, "trials");
config.jets = get_jet_parameters(yaml, "jets");
config.beam = get_Beam(yaml, "beam");
for(size_t i = 0; i < config.process.incoming.size(); ++i){
const auto & in = config.process.incoming[i];
using namespace RHEJ::pid;
if( (in == p || in == p_bar) && in != config.beam.particles[i]){
throw std::invalid_argument{
"Particle type of beam " + std::to_string(i+1) + " incompatible"
+ " with type of incoming particle " + std::to_string(i+1)
};
}
}
set_from_yaml(config.pdf_id, yaml, "pdf");
set_from_yaml(config.unordered_fraction, yaml, "unordered fraction");
if(config.unordered_fraction < 0 || config.unordered_fraction > 1){
throw std::invalid_argument{
"unordered fraction has to be between 0 and 1"
};
}
- config.higgs_properties = get_higgs_properties(yaml, "Higgs properties");
+ config.Higgs_properties = get_higgs_properties(yaml, "Higgs properties");
set_from_yaml_if_defined(config.analysis_parameters, yaml, "analysis");
config.scale_gen = RHEJ::ScaleGenerator{RHEJ::to_ScaleConfig(yaml)};
set_from_yaml_if_defined(config.output, yaml, "event output");
set_from_yaml_if_defined(config.RanLux_init, yaml, "RanLux init");
config.Higgs_coupling = RHEJ::get_Higgs_coupling(yaml, "Higgs coupling");
return config;
}
}
Config load_config(std::string const & config_file){
try{
return to_Config(YAML::LoadFile(config_file));
}
catch(...){
std::cerr << "Error reading " << config_file << ":\n ";
throw;
}
}
}
diff --git a/FixedOrderGen/src/main.cc b/FixedOrderGen/src/main.cc
index 0adaa7d..1b9412c 100644
--- a/FixedOrderGen/src/main.cc
+++ b/FixedOrderGen/src/main.cc
@@ -1,190 +1,191 @@
/**
* Name: main.cc
* Authors: Jeppe R. Andersen
*/
#include <fstream>
#include <algorithm>
#include <memory>
#include <chrono>
#include <iostream>
#include <map>
#include "yaml-cpp/yaml.h"
#include "config.hh"
#include "LHEF/LHEF.h"
#include "RHEJ/CombinedEventWriter.hh"
#include "RHEJ/get_analysis.hh"
#include "RHEJ/utility.hh"
//#include "RHEJ/EventReweighter.hh"
#include "RHEJ/stream.hh"
#include "RHEJ/LesHouchesWriter.hh"
#include "EventGenerator.hh"
#include "PhaseSpacePoint.hh"
namespace{
constexpr auto banner =
" __ ___ __ ______ __ __ \n"
" / / / (_)___ _/ /_ / ____/___ ___ _________ ___ __ / /__ / /______ \n"
" / /_/ / / __ `/ __ \\ / __/ / __ \\/ _ \\/ ___/ __ `/ / / / __ / / _ \\/ __/ ___/ \n"
" / __ / / /_/ / / / / / /___/ / / / __/ / / /_/ / /_/ / / /_/ / __/ /_(__ ) \n"
" /_/ /_/_/\\__, /_/ /_/ /_____/_/ /_/\\___/_/ \\__, /\\__, / \\____/\\___/\\__/____/ \n"
" ____///__/ __ ____ ///__//____/ ______ __ \n"
" / ____(_) _____ ____/ / / __ \\_________/ /__ _____ / ____/__ ____ ___ _________ _/ /_____ _____\n"
" / /_ / / |/_/ _ \\/ __ / / / / / ___/ __ / _ \\/ ___/ / / __/ _ \\/ __ \\/ _ \\/ ___/ __ `/ __/ __ \\/ ___/\n"
" / __/ / /> </ __/ /_/ / / /_/ / / / /_/ / __/ / / /_/ / __/ / / / __/ / / /_/ / /_/ /_/ / / \n"
" /_/ /_/_/|_|\\___/\\__,_/ \\____/_/ \\__,_/\\___/_/ \\____/\\___/_/ /_/\\___/_/ \\__,_/\\__/\\____/_/ \n"
;
constexpr double invGeV2_to_pb = 389379292.;
}
HEJFOG::Config load_config(char const * filename){
try{
return HEJFOG::load_config(filename);
}
catch(std::exception const & exc){
std::cerr << "Error: " << exc.what() << '\n';
std::exit(EXIT_FAILURE);
}
}
std::unique_ptr<RHEJ::Analysis> get_analysis(
YAML::Node const & parameters
){
try{
return RHEJ::get_analysis(parameters);
}
catch(std::exception const & exc){
std::cerr << "Failed to load analysis: " << exc.what() << '\n';
std::exit(EXIT_FAILURE);
}
}
std::string progress_bar(){
std::string result = "0% ";
for(int i = 10; i <= 100; i+= 10){
result += " " + std::to_string(i) + "%";
}
result += "\n|";
for(int i = 10; i <= 100; i+= 10){
result += "---------|";
}
return result + '\n';
}
int main(int argn, char** argv) {
if (argn < 2) {
std::cerr << "\n# Usage:\n.FOgen config_file\n";
return EXIT_FAILURE;
}
std::cout << banner;
fastjet::ClusterSequence::print_banner();
using clock = std::chrono::system_clock;
const auto start_time = clock::now();
// read configuration
auto config = load_config(argv[1]);
std::unique_ptr<RHEJ::Analysis> analysis = get_analysis(
config.analysis_parameters
);
assert(analysis != nullptr);
HEJFOG::EventGenerator generator{
config.process,
config.beam,
config.scale_gen,
config.jets,
config.pdf_id,
config.unordered_fraction,
+ config.Higgs_properties,
config.Higgs_coupling
};
//TODO: fix Les Houches init block (HEPRUP)
LHEF::HEPRUP lhefheprup;
lhefheprup.IDBMUP=std::pair<long,long>(config.beam.particles[0], config.beam.particles[1]);
lhefheprup.EBMUP=std::make_pair(config.beam.energy, config.beam.energy);
lhefheprup.PDFGUP=std::make_pair(0,0);
lhefheprup.PDFSUP=std::make_pair(config.pdf_id,config.pdf_id);
lhefheprup.NPRUP=1;
lhefheprup.XSECUP=std::vector<double>(1.);
lhefheprup.XERRUP=std::vector<double>(1.);
lhefheprup.LPRUP=std::vector<int>{1};
RHEJ::CombinedEventWriter writer{config.output, lhefheprup};
if(config.RanLux_init){
HEJFOG::PhaseSpacePoint::reset_ranlux(*config.RanLux_init);
}
std::map<HEJFOG::Status, int> status_counter;
// Perform N trial generations
int iprint = 0;
double FKL_xs = 0., FKL_xs_err = 0.;
double uno_xs = 0., uno_xs_err = 0.;
std::cout << '\n' << progress_bar();
for (int trials = 0; trials< config.trials;trials++){
if (trials==iprint) {
std::cout << ".";
std::cout.flush();
iprint+=(int) config.trials/100;
}
auto ev = generator.gen_event();
++status_counter[generator.status()];
if(generator.status() != HEJFOG::good) continue;
ev.central().weight /= config.trials;
ev.central().weight *= invGeV2_to_pb;
for(auto & var: ev.variations()){
var.weight /= config.trials;
var.weight *= invGeV2_to_pb;
}
if(analysis->pass_cuts(ev)){
analysis->fill(ev);
writer.write(ev);
if(ev.type() == RHEJ::event_type::FKL){
FKL_xs += ev.central().weight;
FKL_xs_err += ev.central().weight*ev.central().weight;
}
else {
uno_xs += ev.central().weight;
uno_xs_err += ev.central().weight*ev.central().weight;
}
}
} // main event loop
std::cout << std::endl;
std::chrono::duration<double> run_time = (clock::now() - start_time);
std::cout << "\nTask Runtime: " << run_time.count() << " seconds.\n";
std::cout.precision(10);
std::cout.width(12);
std::cout.setf(std::ios::fixed);
std::cout
<< "\nCross section: " << (FKL_xs + uno_xs)
<< " +- " << std::sqrt(FKL_xs_err + uno_xs_err) << " (pb)"
"\n FKL: " << FKL_xs << " +- " << std::sqrt(FKL_xs_err) << " (pb)"
"\n unordered: " << uno_xs << " +- " << std::sqrt(uno_xs_err) << " (pb)"
"\n";
std::cout << '\n';
for(auto && entry: status_counter){
const double fraction = (double) entry.second/config.trials;
const int percent = std::round(100*fraction);
std::cout << "status "
<< std::left << std::setw(16) << (to_string(entry.first) + ":")
<< " [";
for(int i = 0; i < percent/2; ++i) std::cout << '#';
for(int i = percent/2; i < 50; ++i) std::cout << ' ';
std::cout << "] " << percent << "%\n";
}
}
diff --git a/FixedOrderGen/t/h_2j.cc b/FixedOrderGen/t/h_2j.cc
index b2fd856..a61b7eb 100644
--- a/FixedOrderGen/t/h_2j.cc
+++ b/FixedOrderGen/t/h_2j.cc
@@ -1,56 +1,57 @@
#ifdef NDEBUG
#undef NDEBUG
#endif
#include <algorithm>
#include <cmath>
#include <cassert>
#include <iostream>
#include "config.hh"
#include "EventGenerator.hh"
#include "RHEJ/PDF.hh"
#include "RHEJ/MatrixElement.hh"
using namespace HEJFOG;
int main(){
constexpr double invGeV2_to_pb = 389379292.;
constexpr double xs_ref = 2.0304; //calculated with "old" HEJ svn r3364
auto config = load_config("config_h_2j.yml");
HEJFOG::EventGenerator generator{
config.process,
config.beam,
config.scale_gen,
config.jets,
config.pdf_id,
config.unordered_fraction,
+ config.Higgs_properties,
config.Higgs_coupling
};
double xs = 0., xs_err = 0.;
for (int trials = 0; trials < config.trials; ++trials){
auto ev = generator.gen_event();
if(generator.status() != good) continue;
ev.central().weight *= invGeV2_to_pb;
ev.central().weight /= config.trials;
const auto the_Higgs = std::find_if(
begin(ev.outgoing()), end(ev.outgoing()),
[](RHEJ::Sparticle const & p){ return p.type == RHEJ::ParticleID::h; }
);
assert(the_Higgs != end(ev.outgoing()));
if(std::abs(the_Higgs->rapidity()) > 5.) continue;
xs += ev.central().weight;
xs_err += ev.central().weight*ev.central().weight;
}
xs_err = std::sqrt(xs_err);
std::cout << xs_ref << " ~ " << xs << " +- " << xs_err << '\n';
assert(std::abs(xs - xs_ref) < 3*xs_err);
assert(xs_err < 0.01*xs);
}
diff --git a/FixedOrderGen/t/h_3j.cc b/FixedOrderGen/t/h_3j.cc
index b3aebb0..eea79c2 100644
--- a/FixedOrderGen/t/h_3j.cc
+++ b/FixedOrderGen/t/h_3j.cc
@@ -1,57 +1,58 @@
#ifdef NDEBUG
#undef NDEBUG
#endif
#include <algorithm>
#include <cmath>
#include <cassert>
#include <iostream>
#include "config.hh"
#include "EventGenerator.hh"
#include "RHEJ/PDF.hh"
#include "RHEJ/MatrixElement.hh"
using namespace HEJFOG;
int main(){
constexpr double invGeV2_to_pb = 389379292.;
constexpr double xs_ref = 1.0888; //calculated with "old" HEJ svn r3364
auto config = load_config("config_h_2j.yml");
config.process.njets = 3;
HEJFOG::EventGenerator generator{
config.process,
config.beam,
config.scale_gen,
config.jets,
config.pdf_id,
config.unordered_fraction,
+ config.Higgs_properties,
config.Higgs_coupling
};
double xs = 0., xs_err = 0.;
for (int trials = 0; trials < config.trials; ++trials){
auto ev = generator.gen_event();
if(generator.status() != good) continue;
ev.central().weight *= invGeV2_to_pb;
ev.central().weight /= config.trials;
const auto the_Higgs = std::find_if(
begin(ev.outgoing()), end(ev.outgoing()),
[](RHEJ::Sparticle const & p){ return p.type == RHEJ::ParticleID::h; }
);
assert(the_Higgs != end(ev.outgoing()));
if(std::abs(the_Higgs->rapidity()) > 5.) continue;
xs += ev.central().weight;
xs_err += ev.central().weight*ev.central().weight;
}
xs_err = std::sqrt(xs_err);
std::cout << xs_ref << " ~ " << xs << " +- " << xs_err << '\n';
assert(std::abs(xs - xs_ref) < 3*xs_err);
assert(xs_err < 0.02*xs);
}
diff --git a/FixedOrderGen/t/h_3j_uno1.cc b/FixedOrderGen/t/h_3j_uno1.cc
index aa01a12..3ccb544 100644
--- a/FixedOrderGen/t/h_3j_uno1.cc
+++ b/FixedOrderGen/t/h_3j_uno1.cc
@@ -1,60 +1,61 @@
#ifdef NDEBUG
#undef NDEBUG
#endif
// check that adding uno emissions doesn't change the FKL cross section
#include <algorithm>
#include <cmath>
#include <cassert>
#include <iostream>
#include "config.hh"
#include "EventGenerator.hh"
#include "RHEJ/PDF.hh"
#include "RHEJ/MatrixElement.hh"
using namespace HEJFOG;
int main(){
constexpr double invGeV2_to_pb = 389379292.;
constexpr double xs_ref = 0.02603824; //calculated with "old" HEJ svn r3364
auto config = load_config("config_h_2j.yml");
config.process.njets = 3;
config.process.incoming = {RHEJ::pid::u, RHEJ::pid::u};
- config.unordered_fraction = 0.43;
+ config.unordered_fraction = 0.37;
HEJFOG::EventGenerator generator{
config.process,
config.beam,
config.scale_gen,
config.jets,
config.pdf_id,
config.unordered_fraction,
+ config.Higgs_properties,
config.Higgs_coupling
};
double xs = 0., xs_err = 0.;
int uno_found = 0;
for (int trials = 0; trials < config.trials; ++trials){
auto ev = generator.gen_event();
if(generator.status() != good) continue;
if(ev.type() != RHEJ::event_type::FKL){
++uno_found;
continue;
}
ev.central().weight *= invGeV2_to_pb;
ev.central().weight /= config.trials;
xs += ev.central().weight;
xs_err += ev.central().weight*ev.central().weight;
}
xs_err = std::sqrt(xs_err);
std::cout << xs_ref << " ~ " << xs << " +- " << xs_err << '\n';
std::cout << uno_found << " events with unordered emission\n";
assert(uno_found > 0);
assert(std::abs(xs - xs_ref) < 3*xs_err);
assert(xs_err < 0.05*xs);
}
diff --git a/FixedOrderGen/t/h_3j_uno2.cc b/FixedOrderGen/t/h_3j_uno2.cc
index 8794081..09eebd4 100644
--- a/FixedOrderGen/t/h_3j_uno2.cc
+++ b/FixedOrderGen/t/h_3j_uno2.cc
@@ -1,54 +1,55 @@
#ifdef NDEBUG
#undef NDEBUG
#endif
// check uno cross section
#include <algorithm>
#include <cmath>
#include <cassert>
#include <iostream>
#include "config.hh"
#include "EventGenerator.hh"
#include "RHEJ/PDF.hh"
#include "RHEJ/MatrixElement.hh"
using namespace HEJFOG;
int main(){
constexpr double invGeV2_to_pb = 389379292.;
constexpr double xs_ref = 0.00336012;
auto config = load_config("config_h_2j.yml");
config.process.njets = 3;
config.process.incoming = {RHEJ::pid::u, RHEJ::pid::u};
config.unordered_fraction = 1.;
HEJFOG::EventGenerator generator{
config.process,
config.beam,
config.scale_gen,
config.jets,
config.pdf_id,
config.unordered_fraction,
+ config.Higgs_properties,
config.Higgs_coupling
};
double xs = 0., xs_err = 0.;
for (int trials = 0; trials < config.trials; ++trials){
auto ev = generator.gen_event();
if(ev.type() == RHEJ::event_type::FKL) continue;
if(generator.status() != good) continue;
ev.central().weight *= invGeV2_to_pb;
ev.central().weight /= config.trials;
xs += ev.central().weight;
xs_err += ev.central().weight*ev.central().weight;
}
xs_err = std::sqrt(xs_err);
std::cout << xs_ref << " ~ " << xs << " +- " << xs_err << '\n';
assert(std::abs(xs - xs_ref) < 3*xs_err);
assert(xs_err < 0.05*xs);
}
diff --git a/FixedOrderGen/t/h_5j.cc b/FixedOrderGen/t/h_5j.cc
index d8e4f6f..1cdd2fa 100644
--- a/FixedOrderGen/t/h_5j.cc
+++ b/FixedOrderGen/t/h_5j.cc
@@ -1,52 +1,53 @@
#ifdef NDEBUG
#undef NDEBUG
#endif
// This is a regression test
// the reference cross section has not been checked against any other program
#include <algorithm>
#include <cmath>
#include <cassert>
#include <iostream>
#include "config.hh"
#include "EventGenerator.hh"
#include "RHEJ/PDF.hh"
#include "RHEJ/MatrixElement.hh"
using namespace HEJFOG;
int main(){
constexpr double invGeV2_to_pb = 389379292.;
constexpr double xs_ref = 0.235969;
auto config = load_config("config_h_2j.yml");
config.process.njets = 5;
HEJFOG::EventGenerator generator{
config.process,
config.beam,
config.scale_gen,
config.jets,
config.pdf_id,
config.unordered_fraction,
+ config.Higgs_properties,
config.Higgs_coupling
};
double xs = 0., xs_err = 0.;
for (int trials = 0; trials < config.trials; ++trials){
auto ev = generator.gen_event();
if(generator.status() != good) continue;
ev.central().weight *= invGeV2_to_pb;
ev.central().weight /= config.trials;
xs += ev.central().weight;
xs_err += ev.central().weight*ev.central().weight;
}
xs_err = std::sqrt(xs_err);
std::cout << xs_ref << " ~ " << xs << " +- " << xs_err << '\n';
assert(std::abs(xs - xs_ref) < 3*xs_err);
assert(xs_err < 0.06*xs);
}
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Tue, Nov 19, 6:10 PM (1 d, 19 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3804283
Default Alt Text
(48 KB)
Attached To
rHEJ HEJ
Event Timeline
Log In to Comment