Page MenuHomeHEPForge

No OneTemporary

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

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)

Event Timeline