Page MenuHomeHEPForge
Files F9501704

No OneTemporary
Actions

View Options

diff --git a/Changes-API.md b/Changes-API.md
index d96807a..60baa9e 100644
--- a/Changes-API.md
+++ b/Changes-API.md
@@ -1,94 +1,99 @@
# Changelog for HEJ API
This log lists only changes on the HEJ API. These are primarily code changes
relevant for calling HEJ as an API. This file should only be read as an addition
to [`Changes.md`](Changes.md), where the main features are documented.
## Version 2.1
### 2.1.0
#### Changes to Event class
* Restructured `Event` class
- `Event` can now only be build from a (new) `Event::EventData` class
- Removed default constructor for `Event`
- `Event::EventData` replaces the old `UnclusteredEvent` struct.
- `UnclusteredEvent` is now deprecated, and will be removed in HEJ Version
2.3.0
- Removed `Event.unclustered()` function
- Added new member function `Events.parameters()`, to directly access
(underlying) `Parameters<EventParameters>`
- New member functions `begin_partons`, `end_partons` (`rbegin_partons`,
`rend_partons`) with aliases `cbegin_partons`, `cend_partons`
(`crbegin_partons`, `crend_partons`) for constant (reversed) iterators over
outgoing partons.
* New function `Event::EventData.reconstruct_intermediate()` to reconstruct
bosons from decays, e.g. `positron + nu_e => Wp`
* Added optional Colour charges to particles (`Particle.colour`)
- Colour connection in the HEJ limit can be generated via
`Event.generate_colours` (automatically done in the resummation)
- Colour configuration of input events can be checked with
`Event.is_leading_colour`
* Added function `Event.valid_hej_state` to check if event _could have_ been
produced by `HEJ` according to the `max ext soft pt fraction` cut on the jets
#### Updated functions
* Renamed `EventType::nonHEJ` to `EventType::non_resummable` and `is_HEJ()`
to `is_resummable()` such that run card is consistent with internal workings
* Made `MatrixElement.tree_kin(...)` and `MatrixElement.tree_param(...)` public
* New `EventReweighter` member function `treatment` to query the
treatment with respect to resummation for the various event types.
* Added auxiliary functions to obtain a `std::string` from and
`EventDescription` (`to_string` for human readable, and `to_simple_string`
for easy parsable string)
#### New classes
* New class `Unweighter` to do unweighting
* New class `CrossSectionAccumulator` to keep track of Cross Section of the
different subprocess
* New template struct `Parameters` similar to old `Weights`
- `Weights` are now an alias for `Parameters<double>`. Calling `Weights` did
not change
- `Weights.hh` was replaced by `Parameters.hh`. The old `Weights.hh` header
will be removed in HEJ Version 2.2.0
* Function to multiplication and division of `EventParameters.weight` by double
- This can be combined with `Parameters`, e.g.
`Parameters<EventParameters>*Weights`, see also `Events.parameters()`
- Moved `EventParameters` to `Parameters.hh` header
* new class `EWConstants` replaces previously hard coded `vev`
- `EWConstants` have to be set in the general `Config` and the
`MatrixElementConfig`
#### Input/Output
* New abstract `EventReader` class, as base for reading events from files
- Moved LHE file reader to `HEJ::LesHouchesReader`
* Support reading (`HDF5Reader`) and writing (`HDF5Writer`) `hdf5` files
* New `BufferedEventReader` class that allows to put events back into
the reader.
* New `SherpaLHEReader` to read Sherpa LHE files with correct weights
* `get_analysis` now requires `YAML::Node` and `LHEF::HEPRUP` as arguments
* Replaced `HepMCInterface` and `HepMCWriter` by `HepMCInterfaceX` and
`HepMCWriterX` respectively, with `X` being the major version of HepMC (2 or
3)
- Renamed `HepMCInterfaceX::init_kinematics` to `HepMCInterfaceX::init_event`
and protected it, use `HepMCInterfaceX::operator()` instead
- Removed redundant `HepMCInterfaceX::add_variation` function
#### Linking
* Export cmake target `HEJ::HEJ` to link directly against `libHEJ`
#### Miscellaneous
* Capitalisation of `Config.hh` now matches class `Config` (was `config.hh`)
* Replaced redundant member `EventReweighterConfig::jet_param` getter function
`EventReweighter.jet_param()` (`JetParameters` are already in
`PhaseSpacePointConfig`)
+* Avoid storing reference to the Random Number Generator inside classes
+ - Constructors of `EventReweighter` now expect `std::shared_ptr<HEJ::RNG>`
+ (was reference)
+ - Moved reference to `HEJ::RNG` from constructor of `JetSplitter` to
+ `JetSplitter.split`
## Version 2.0
### 2.0.4
* Fixed wrong path of `HEJ_INCLUDE_DIR` in `hej-config.cmake`
### 2.0.0
* First release
diff --git a/FixedOrderGen/include/EventGenerator.hh b/FixedOrderGen/include/EventGenerator.hh
index e345a81..b89d640 100644
--- a/FixedOrderGen/include/EventGenerator.hh
+++ b/FixedOrderGen/include/EventGenerator.hh
@@ -1,65 +1,67 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
+#include <memory>
+
#include "HEJ/MatrixElement.hh"
#include "HEJ/optional.hh"
#include "HEJ/PDF.hh"
#include "HEJ/RNG.hh"
#include "Beam.hh"
#include "Decay.hh"
#include "JetParameters.hh"
#include "Process.hh"
#include "Status.hh"
namespace HEJ{
class Event;
class HiggsCouplingSettings;
class ScaleGenerator;
class EWConstants;
}
//! Namespace for HEJ Fixed Order Generator
namespace HEJFOG{
class EventGenerator{
public:
EventGenerator(
Process process,
Beam beam,
HEJ::ScaleGenerator scale_gen,
JetParameters jets,
int pdf_id,
double subl_change,
unsigned int subl_channels,
ParticlesDecayMap particle_decays,
HEJ::HiggsCouplingSettings Higgs_coupling,
HEJ::EWConstants ew_parameters,
- HEJ::RNG & ran
+ std::shared_ptr<HEJ::RNG> ran
);
HEJ::optional<HEJ::Event> gen_event();
Status status() const {
return status_;
}
private:
HEJ::PDF pdf_;
HEJ::MatrixElement ME_;
HEJ::ScaleGenerator scale_gen_;
Process process_;
JetParameters jets_;
Beam beam_;
Status status_;
double subl_change_;
unsigned int subl_channels_;
ParticlesDecayMap particle_decays_;
HEJ::EWConstants ew_parameters_;
- std::reference_wrapper<HEJ::RNG> ran_;
+ std::shared_ptr<HEJ::RNG> ran_;
};
}
diff --git a/FixedOrderGen/src/EventGenerator.cc b/FixedOrderGen/src/EventGenerator.cc
index c92661e..58a5d62 100644
--- a/FixedOrderGen/src/EventGenerator.cc
+++ b/FixedOrderGen/src/EventGenerator.cc
@@ -1,93 +1,95 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "EventGenerator.hh"
+#include <utility>
+
#include "Process.hh"
#include "Beam.hh"
#include "JetParameters.hh"
#include "PhaseSpacePoint.hh"
#include "HEJ/Config.hh"
#include "HEJ/Event.hh"
#include "HEJ/EWConstants.hh"
namespace HEJFOG{
EventGenerator::EventGenerator(
Process process,
Beam beam,
HEJ::ScaleGenerator scale_gen,
JetParameters jets,
int pdf_id,
double subl_change,
unsigned int subl_channels,
ParticlesDecayMap particle_decays,
HEJ::HiggsCouplingSettings Higgs_coupling,
HEJ::EWConstants ew_parameters,
- HEJ::RNG & ran
+ std::shared_ptr<HEJ::RNG> ran
):
pdf_{pdf_id, beam.particles[0], beam.particles[1]},
ME_{
[this](double mu){ return pdf_.Halphas(mu); },
HEJ::MatrixElementConfig{
false,
std::move(Higgs_coupling),
ew_parameters
}
},
scale_gen_{std::move(scale_gen)},
process_{std::move(process)},
jets_{std::move(jets)},
beam_{std::move(beam)},
subl_change_{subl_change},
subl_channels_{subl_channels},
particle_decays_{std::move(particle_decays)},
ew_parameters_{ew_parameters},
- ran_{ran}
+ ran_{std::move(ran)}
{
}
HEJ::optional<HEJ::Event> EventGenerator::gen_event(){
HEJFOG::PhaseSpacePoint psp{
process_,
jets_,
pdf_, beam_.energy,
subl_change_, subl_channels_,
particle_decays_,
ew_parameters_,
- ran_
+ *ran_
};
status_ = psp.status();
if(status_ != good) return {};
HEJ::Event ev = scale_gen_(
HEJ::Event{
to_EventData( std::move(psp) ).cluster( jets_.def, jets_.min_pt)
}
);
if(!is_resummable(ev.type()))
throw HEJ::not_implemented("Tried to generate a event type, "
"which is not yet implemented in HEJ.");
- ev.generate_colours(ran_);
+ ev.generate_colours(*ran_);
const double shat = HEJ::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
ev.parameters() *= ME_.tree(ev)/(shat*shat);
// and PDFs
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(size_t i = 0; i < ev.variations().size(); ++i){
auto & var = ev.variations(i);
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/main.cc b/FixedOrderGen/src/main.cc
index 22eaa3f..3451bdc 100644
--- a/FixedOrderGen/src/main.cc
+++ b/FixedOrderGen/src/main.cc
@@ -1,260 +1,261 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include <algorithm>
#include <chrono>
#include <fstream>
#include <iostream>
#include <map>
#include <memory>
#include <cstdint>
#include "LHEF/LHEF.h"
#include "yaml-cpp/yaml.h"
#include <boost/iterator/filter_iterator.hpp>
#include "HEJ/CombinedEventWriter.hh"
#include "HEJ/CrossSectionAccumulator.hh"
#include "HEJ/get_analysis.hh"
#include "HEJ/LesHouchesWriter.hh"
#include "HEJ/make_RNG.hh"
#include "HEJ/ProgressBar.hh"
#include "HEJ/stream.hh"
#include "HEJ/Unweighter.hh"
#include "config.hh"
#include "EventGenerator.hh"
#include "PhaseSpacePoint.hh"
#include "Version.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<HEJ::Analysis> get_analysis(
YAML::Node const & parameters, LHEF::HEPRUP const & heprup
){
try{
return HEJ::get_analysis(parameters, heprup);
}
catch(std::exception const & exc){
std::cerr << "Failed to load analysis: " << exc.what() << '\n';
std::exit(EXIT_FAILURE);
}
}
template<class Iterator>
auto make_lowpt_filter(Iterator begin, Iterator end, HEJ::optional<double> peak_pt){
return boost::make_filter_iterator(
[peak_pt](HEJ::Event const & ev){
assert(! ev.jets().empty());
double min_pt = peak_pt?(*peak_pt):0.;
const auto softest_jet = fastjet::sorted_by_pt(ev.jets()).back();
return softest_jet.pt() > min_pt;
},
begin, end
);
}
int main(int argn, char** argv) {
using namespace std::string_literals;
if (argn < 2) {
std::cerr << "\n# Usage:\n." << argv[0] << " config_file\n";
return EXIT_FAILURE;
}
std::cout << banner;
std::cout << "Version " << HEJFOG::Version::String()
<< ", revision " << HEJFOG::Version::revision() << std::endl;
fastjet::ClusterSequence::print_banner();
using clock = std::chrono::system_clock;
const auto start_time = clock::now();
// read configuration
auto config = load_config(argv[1]);
- auto ran = HEJ::make_RNG(config.rng.name, config.rng.seed);
+ std::shared_ptr<HEJ::RNG> ran{
+ HEJ::make_RNG(config.rng.name, config.rng.seed)};
assert(ran != nullptr);
HEJ::ScaleGenerator scale_gen{
config.scales.base,
config.scales.factors,
config.scales.max_ratio
};
HEJFOG::EventGenerator generator{
config.process,
config.beam,
std::move(scale_gen),
config.jets,
config.pdf_id,
config.subleading_fraction,
config.subleading_channels,
config.particle_decays,
config.Higgs_coupling,
config.ew_parameters,
- *ran
+ ran
};
// prepare process information for output
LHEF::HEPRUP heprup;
heprup.IDBMUP=std::pair<long,long>(config.beam.particles[0], config.beam.particles[1]);
heprup.EBMUP=std::make_pair(config.beam.energy, config.beam.energy);
heprup.PDFGUP=std::make_pair(0,0);
heprup.PDFSUP=std::make_pair(config.pdf_id,config.pdf_id);
heprup.NPRUP=1;
heprup.XSECUP=std::vector<double>(1.);
heprup.XERRUP=std::vector<double>(1.);
heprup.LPRUP=std::vector<int>{1};
heprup.generators.emplace_back(LHEF::XMLTag{});
heprup.generators.back().name = HEJFOG::Version::package_name();
heprup.generators.back().version = HEJFOG::Version::String();
HEJ::CombinedEventWriter writer{config.output, heprup};
std::unique_ptr<HEJ::Analysis> analysis = get_analysis(
config.analysis_parameters, heprup
);
assert(analysis != nullptr);
// warm-up phase to train unweighter
HEJ::optional<HEJ::Unweighter> unweighter{};
std::map<HEJFOG::Status, std::uint64_t> status_counter;
std::vector<HEJ::Event> events;
std::uint64_t trials = 0;
if(config.unweight) {
std::cout << "Calibrating unweighting ...\n";
const auto warmup_start = clock::now();
const size_t warmup_events = config.unweight->sample_size;
HEJ::ProgressBar<size_t> warmup_progress{std::cout, warmup_events};
for(; events.size() < warmup_events; ++trials){
auto ev = generator.gen_event();
++status_counter[generator.status()];
assert( (generator.status() == HEJFOG::good) == bool(ev) );
if(generator.status() == HEJFOG::good && analysis->pass_cuts(*ev, *ev)) {
events.emplace_back(std::move(*ev));
++warmup_progress;
}
}
std::cout << std::endl;
unweighter = HEJ::Unweighter();
unweighter->set_cut_to_peakwt(
make_lowpt_filter(events.cbegin(), events.cend(), config.jets.peak_pt),
make_lowpt_filter(events.cend(), events.cend(), config.jets.peak_pt),
config.unweight->max_dev
);
std::vector<HEJ::Event> unweighted_events;
for(auto && ev: events) {
auto unweighted = unweighter->unweight(std::move(ev), *ran);
if(unweighted) {
unweighted_events.emplace_back(std::move(*unweighted));
}
}
events = std::move(unweighted_events);
if(events.empty()) {
std::cerr <<
"Failed to generate events. Please increase \"unweight: sample size\""
" or reduce \"unweight: max deviation\"\n";
return EXIT_FAILURE;
}
const auto warmup_end = clock::now();
const double completion = static_cast<double>(events.size())/config.events;
const std::chrono::duration<double> remaining_time =
(warmup_end- warmup_start)*(1./completion - 1);
const auto finish = clock::to_time_t(
std::chrono::time_point_cast<std::chrono::seconds>(warmup_end + remaining_time)
);
std::cout
<< "Generated " << events.size() << "/" << config.events << " events ("
<< static_cast<int>(std::round(100*completion)) << "%)\n"
<< "Estimated remaining generation time: "
<< remaining_time.count() << " seconds ("
<< std::put_time(std::localtime(&finish), "%c") << ")\n\n";
} // end unweighting warm-up
// main generation loop
// event weight is wrong, need to divide by "total number of trials" afterwards
HEJ::ProgressBar<size_t> progress{std::cout, config.events};
progress.increment(events.size());
events.reserve(config.events);
for(; events.size() < config.events; ++trials){
auto ev = generator.gen_event();
++status_counter[generator.status()];
assert( (generator.status() == HEJFOG::good) == bool(ev) );
if(generator.status() == HEJFOG::good && analysis->pass_cuts(*ev, *ev)) {
if(unweighter) {
auto unweighted = unweighter->unweight(std::move(*ev), *ran);
if(! unweighted) continue;
ev = std::move(unweighted);
}
events.emplace_back(std::move(*ev));
++progress;
}
}
std::cout << std::endl;
// final run though events with correct weight
HEJ::CrossSectionAccumulator xs;
for(auto & ev: events){
ev.parameters() *= invGeV2_to_pb/trials;
analysis->fill(ev, ev);
writer.write(ev);
xs.fill(ev);
}
analysis->finalise();
// Print final informations
const std::chrono::duration<double> run_time = (clock::now() - start_time);
std::cout << "\nTask Runtime: " << run_time.count() << " seconds for "
<< events.size() << " Events (" << events.size()/run_time.count()
<< " evts/s)\n" << std::endl;
std::cout << xs << "\n";
for(auto && entry: status_counter){
const double fraction = static_cast<double>(entry.second)/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 << "%" << std::endl;
}
return EXIT_SUCCESS;
}
diff --git a/FixedOrderGen/t/2j.cc b/FixedOrderGen/t/2j.cc
index 63419e9..498e352 100644
--- a/FixedOrderGen/t/2j.cc
+++ b/FixedOrderGen/t/2j.cc
@@ -1,67 +1,67 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#ifdef NDEBUG
#undef NDEBUG
#endif
#include <algorithm>
#include <cmath>
#include <cassert>
#include <iostream>
#include "config.hh"
#include "EventGenerator.hh"
#include "HEJ/Mixmax.hh"
#include "HEJ/Event.hh"
#include "HEJ/PDF.hh"
#include "HEJ/MatrixElement.hh"
using namespace HEJFOG;
int main(){
constexpr double invGeV2_to_pb = 389379292.;
constexpr double xs_ref = 86.42031848*1e6; //calculated with "combined" HEJ svn r3480
auto config = load_config("config_2j.yml");
- HEJ::Mixmax ran{};
+ std::shared_ptr<HEJ::RNG> ran{std::make_shared<HEJ::Mixmax>()};
HEJFOG::EventGenerator generator{
config.process,
config.beam,
HEJ::ScaleGenerator{
config.scales.base,
config.scales.factors,
config.scales.max_ratio
},
config.jets,
config.pdf_id,
config.subleading_fraction,
config.subleading_channels,
config.particle_decays,
config.Higgs_coupling,
config.ew_parameters,
ran
};
double xs = 0., xs_err = 0.;
for (size_t trials = 0; trials < config.events; ++trials){
auto ev = generator.gen_event();
if(generator.status() != good) continue;
assert(ev);
ev->central().weight *= invGeV2_to_pb;
ev->central().weight /= config.events;
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 << std::endl;
assert(std::abs(xs - xs_ref) < 3*xs_err);
assert(xs_err < 0.01*xs);
return EXIT_SUCCESS;
}
diff --git a/FixedOrderGen/t/4j.cc b/FixedOrderGen/t/4j.cc
index e3432e5..c2b4634 100644
--- a/FixedOrderGen/t/4j.cc
+++ b/FixedOrderGen/t/4j.cc
@@ -1,69 +1,69 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#ifdef NDEBUG
#undef NDEBUG
#endif
#include <algorithm>
#include <cmath>
#include <cassert>
#include <iostream>
#include "config.hh"
#include "EventGenerator.hh"
#include "HEJ/Mixmax.hh"
#include "HEJ/Event.hh"
#include "HEJ/PDF.hh"
#include "HEJ/MatrixElement.hh"
using namespace HEJFOG;
int main(){
constexpr double invGeV2_to_pb = 389379292.;
// calculated with 13207b5f67a5f40a2141aa7ee515b022bd4efb65
constexpr double xs_ref = 915072; //+- 7227.72
auto config = load_config("config_2j.yml");
config.process.njets = 4;
- HEJ::Mixmax ran{};
+ std::shared_ptr<HEJ::RNG> ran{std::make_shared<HEJ::Mixmax>()};
HEJFOG::EventGenerator generator{
config.process,
config.beam,
HEJ::ScaleGenerator{
config.scales.base,
config.scales.factors,
config.scales.max_ratio
},
config.jets,
config.pdf_id,
config.subleading_fraction,
config.subleading_channels,
config.particle_decays,
config.Higgs_coupling,
config.ew_parameters,
ran
};
double xs = 0., xs_err = 0.;
for (size_t trials = 0; trials < config.events; ++trials){
auto ev = generator.gen_event();
if(generator.status() != good) continue;
assert(ev);
ev->central().weight *= invGeV2_to_pb;
ev->central().weight /= config.events;
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 << std::endl;
assert(std::abs(xs - xs_ref) < 3*xs_err);
assert(xs_err < 0.03*xs);
return EXIT_SUCCESS;
}
diff --git a/FixedOrderGen/t/W_reconstruct_enu.cc b/FixedOrderGen/t/W_reconstruct_enu.cc
index d3918c8..150e3c3 100644
--- a/FixedOrderGen/t/W_reconstruct_enu.cc
+++ b/FixedOrderGen/t/W_reconstruct_enu.cc
@@ -1,72 +1,72 @@
/**
* \brief that the reconstruction of the W works
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#ifdef NDEBUG
#undef NDEBUG
#endif
#include <algorithm>
#include "config.hh"
#include "EventGenerator.hh"
#include "HEJ/Event.hh"
#include "HEJ/Mixmax.hh"
using namespace HEJFOG;
using namespace HEJ;
namespace {
constexpr size_t num_events = 1000;
constexpr double invGeV2_to_pb = 389379292.;
}
double get_xs(std::string config_name){
auto config { load_config(config_name) };
config.events = num_events;
- HEJ::Mixmax ran{1};
+ std::shared_ptr<HEJ::RNG> ran{std::make_shared<HEJ::Mixmax>(11)};
HEJFOG::EventGenerator generator{
config.process,
config.beam,
HEJ::ScaleGenerator{
config.scales.base,
config.scales.factors,
config.scales.max_ratio
},
config.jets,
config.pdf_id,
config.subleading_fraction,
config.subleading_channels,
config.particle_decays,
config.Higgs_coupling,
config.ew_parameters,
ran
};
double xs = 0.;
for (size_t trials = 0; trials < config.events; ++trials){
auto ev = generator.gen_event();
if(generator.status() != good) continue;
assert(ev);
ev->central().weight *= invGeV2_to_pb;
ev->central().weight /= config.events;
xs += ev->central().weight;
}
return xs;
}
int main(){
double xs_W{ get_xs("config_Wp_2j.yml")};
double xs_enu{ get_xs("config_Wp_2j_decay.yml")};
if(std::abs(xs_W/xs_enu-1.)>1e-6){
std::cerr << "Reconstructing the W in the runcard gave a different results ("
<< xs_W << " vs. "<< xs_enu << " -> " << std::abs(xs_W/xs_enu-1.)*100 << "%)\n";
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
diff --git a/FixedOrderGen/t/h_2j.cc b/FixedOrderGen/t/h_2j.cc
index c881ce9..84829b8 100644
--- a/FixedOrderGen/t/h_2j.cc
+++ b/FixedOrderGen/t/h_2j.cc
@@ -1,75 +1,75 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#ifdef NDEBUG
#undef NDEBUG
#endif
#include <algorithm>
#include <cmath>
#include <cassert>
#include <iostream>
#include "config.hh"
#include "EventGenerator.hh"
#include "HEJ/Mixmax.hh"
#include "HEJ/Event.hh"
#include "HEJ/PDF.hh"
#include "HEJ/MatrixElement.hh"
using namespace HEJFOG;
int main(){
constexpr double invGeV2_to_pb = 389379292.;
constexpr double xs_ref = 2.04928; // +- 0.00377252
//calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20
auto config = load_config("config_h_2j.yml");
- HEJ::Mixmax ran{};
+ std::shared_ptr<HEJ::RNG> ran{std::make_shared<HEJ::Mixmax>()};
HEJFOG::EventGenerator generator{
config.process,
config.beam,
HEJ::ScaleGenerator{
config.scales.base,
config.scales.factors,
config.scales.max_ratio
},
config.jets,
config.pdf_id,
config.subleading_fraction,
config.subleading_channels,
config.particle_decays,
config.Higgs_coupling,
config.ew_parameters,
ran
};
double xs = 0., xs_err = 0.;
for (size_t trials = 0; trials < config.events; ++trials){
auto ev = generator.gen_event();
if(generator.status() != good) continue;
assert(ev);
ev->central().weight *= invGeV2_to_pb;
ev->central().weight /= config.events;
const auto the_Higgs = std::find_if(
begin(ev->outgoing()), end(ev->outgoing()),
[](HEJ::Particle const & p){ return p.type == HEJ::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 << std::endl;
assert(std::abs(xs - xs_ref) < 3*xs_err);
assert(xs_err < 0.01*xs);
return EXIT_SUCCESS;
}
diff --git a/FixedOrderGen/t/h_2j_decay.cc b/FixedOrderGen/t/h_2j_decay.cc
index 3dcac48..163e74d 100644
--- a/FixedOrderGen/t/h_2j_decay.cc
+++ b/FixedOrderGen/t/h_2j_decay.cc
@@ -1,94 +1,94 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#ifdef NDEBUG
#undef NDEBUG
#endif
#include <algorithm>
#include <cmath>
#include <cassert>
#include <iostream>
#include "config.hh"
#include "EventGenerator.hh"
#include "HEJ/Event.hh"
#include "HEJ/MatrixElement.hh"
#include "HEJ/Particle.hh"
#include "HEJ/PDF.hh"
#include "HEJ/Ranlux64.hh"
#include "HEJ/utility.hh"
using namespace HEJFOG;
bool pass_dR_cut(
std::vector<fastjet::PseudoJet> const & jets,
std::vector<HEJ::Particle> const & photons
){
constexpr double delta_R_min = 0.7;
for(auto const & jet: jets){
for(auto const & photon: photons){
if(jet.delta_R(photon.p) < delta_R_min) return false;
}
}
return true;
}
int main(){
constexpr double invGeV2_to_pb = 389379292.;
constexpr double xs_ref = 0.00429198; // +- 1.0488e-05
//calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20
auto config = load_config("config_h_2j_decay.yml");
- HEJ::Ranlux64 ran{};
+ std::shared_ptr<HEJ::RNG> ran{std::make_shared<HEJ::Ranlux64>()};
HEJFOG::EventGenerator generator{
config.process,
config.beam,
HEJ::ScaleGenerator{
config.scales.base,
config.scales.factors,
config.scales.max_ratio
},
config.jets,
config.pdf_id,
config.subleading_fraction,
config.subleading_channels,
config.particle_decays,
config.Higgs_coupling,
config.ew_parameters,
ran
};
double xs = 0., xs_err = 0.;
for (size_t trials = 0; trials < config.events; ++trials){
auto ev = generator.gen_event();
if(generator.status() != good) continue;
assert(ev);
assert(ev->decays().size() == 1);
const auto decay = begin(ev->decays());
assert(ev->outgoing().size() > decay->first);
const auto & the_Higgs = ev->outgoing()[decay->first];
assert(the_Higgs.type == HEJ::pid::Higgs);
assert(decay->second.size() == 2);
auto const & gamma = decay->second;
assert(gamma[0].type == HEJ::pid::photon);
assert(gamma[1].type == HEJ::pid::photon);
assert(HEJ::nearby_ep(gamma[0].p + gamma[1].p, the_Higgs.p, 1e-6));
if(!pass_dR_cut(ev->jets(), gamma)) continue;
ev->central().weight *= invGeV2_to_pb;
ev->central().weight /= config.events;
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 << std::endl;
assert(std::abs(xs - xs_ref) < 3*xs_err);
assert(xs_err < 0.012*xs);
return EXIT_SUCCESS;
}
diff --git a/FixedOrderGen/t/h_3j.cc b/FixedOrderGen/t/h_3j.cc
index 292cc62..16a6670 100644
--- a/FixedOrderGen/t/h_3j.cc
+++ b/FixedOrderGen/t/h_3j.cc
@@ -1,76 +1,76 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#ifdef NDEBUG
#undef NDEBUG
#endif
#include <algorithm>
#include <cmath>
#include <cassert>
#include <iostream>
#include "config.hh"
#include "EventGenerator.hh"
#include "HEJ/Ranlux64.hh"
#include "HEJ/Event.hh"
#include "HEJ/MatrixElement.hh"
#include "HEJ/PDF.hh"
using namespace HEJFOG;
int main(){
constexpr double invGeV2_to_pb = 389379292.;
constexpr double xs_ref = 1.07807; // +- 0.0071
//calculated with HEJ revision 93efdc851b02a907a6fcc63956387f9f4c1111c2 +1
auto config = load_config("config_h_2j.yml");
config.process.njets = 3;
- HEJ::Ranlux64 ran{};
+ std::shared_ptr<HEJ::RNG> ran{std::make_shared<HEJ::Ranlux64>()};
HEJFOG::EventGenerator generator{
config.process,
config.beam,
HEJ::ScaleGenerator{
config.scales.base,
config.scales.factors,
config.scales.max_ratio
},
config.jets,
config.pdf_id,
config.subleading_fraction,
config.subleading_channels,
config.particle_decays,
config.Higgs_coupling,
config.ew_parameters,
ran
};
double xs = 0., xs_err = 0.;
for (size_t trials = 0; trials < config.events; ++trials){
auto ev = generator.gen_event();
if(generator.status() != good) continue;
assert(ev);
ev->central().weight *= invGeV2_to_pb;
ev->central().weight /= config.events;
const auto the_Higgs = std::find_if(
begin(ev->outgoing()), end(ev->outgoing()),
[](HEJ::Particle const & p){ return p.type == HEJ::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 << std::endl;
assert(std::abs(xs - xs_ref) < 3*xs_err);
assert(xs_err < 0.02*xs);
return EXIT_SUCCESS;
}
diff --git a/FixedOrderGen/t/h_3j_uno1.cc b/FixedOrderGen/t/h_3j_uno1.cc
index 7e8042b..09fd6eb 100644
--- a/FixedOrderGen/t/h_3j_uno1.cc
+++ b/FixedOrderGen/t/h_3j_uno1.cc
@@ -1,80 +1,80 @@
/**
* check that adding uno emissions doesn't change the FKL cross section
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#ifdef NDEBUG
#undef NDEBUG
#endif
#include <algorithm>
#include <cassert>
#include <cmath>
#include <iostream>
#include "config.hh"
#include "EventGenerator.hh"
#include "HEJ/Ranlux64.hh"
#include "Subleading.hh"
#include "HEJ/Event.hh"
#include "HEJ/MatrixElement.hh"
#include "HEJ/PDF.hh"
using namespace HEJFOG;
int main(){
constexpr double invGeV2_to_pb = 389379292.;
constexpr double xs_ref = 0.0243548; // +- 0.000119862
//calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20
auto config = load_config("config_h_2j.yml");
config.process.njets = 3;
config.process.incoming = {HEJ::pid::u, HEJ::pid::u};
config.subleading_channels = HEJFOG::Subleading::uno;
- HEJ::Ranlux64 ran{};
+ std::shared_ptr<HEJ::RNG> ran{std::make_shared<HEJ::Ranlux64>()};
HEJFOG::EventGenerator generator{
config.process,
config.beam,
HEJ::ScaleGenerator{
config.scales.base,
config.scales.factors,
config.scales.max_ratio
},
config.jets,
config.pdf_id,
config.subleading_fraction,
config.subleading_channels,
config.particle_decays,
config.Higgs_coupling,
config.ew_parameters,
ran
};
double xs = 0., xs_err = 0.;
size_t uno_found = 0;
for (size_t trials = 0; trials < config.events; ++trials){
auto ev = generator.gen_event();
if(generator.status() != good) continue;
assert(ev);
if(ev->type() != HEJ::event_type::FKL){
++uno_found;
continue;
}
ev->central().weight *= invGeV2_to_pb;
ev->central().weight /= config.events;
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" << std::endl;
assert(uno_found > 0);
assert(std::abs(xs - xs_ref) < 3*xs_err);
assert(xs_err < 0.05*xs);
return EXIT_SUCCESS;
}
diff --git a/FixedOrderGen/t/h_3j_uno2.cc b/FixedOrderGen/t/h_3j_uno2.cc
index 3f1c5f8..b93f1ad 100644
--- a/FixedOrderGen/t/h_3j_uno2.cc
+++ b/FixedOrderGen/t/h_3j_uno2.cc
@@ -1,75 +1,75 @@
/**
* check uno cross section
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#ifdef NDEBUG
#undef NDEBUG
#endif
#include <algorithm>
#include <cassert>
#include <cmath>
#include <iostream>
#include "config.hh"
#include "EventGenerator.hh"
#include "HEJ/Ranlux64.hh"
#include "Subleading.hh"
#include "HEJ/Event.hh"
#include "HEJ/MatrixElement.hh"
#include "HEJ/PDF.hh"
using namespace HEJFOG;
int main(){
constexpr double invGeV2_to_pb = 389379292.;
constexpr double xs_ref = 0.00347538; // +- 3.85875e-05
//calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20
auto config = load_config("config_h_2j.yml");
config.process.njets = 3;
config.process.incoming = {HEJ::pid::u, HEJ::pid::u};
config.subleading_fraction = 1.;
config.subleading_channels = HEJFOG::Subleading::uno;
- HEJ::Ranlux64 ran{};
+ std::shared_ptr<HEJ::RNG> ran{std::make_shared<HEJ::Ranlux64>()};
HEJFOG::EventGenerator generator{
config.process,
config.beam,
HEJ::ScaleGenerator{
config.scales.base,
config.scales.factors,
config.scales.max_ratio
},
config.jets,
config.pdf_id,
config.subleading_fraction,
config.subleading_channels,
config.particle_decays,
config.Higgs_coupling,
config.ew_parameters,
ran
};
double xs = 0., xs_err = 0.;
for (size_t trials = 0; trials < config.events; ++trials){
auto ev = generator.gen_event();
if(generator.status() != good) continue;
assert(ev);
if(ev->type() == HEJ::event_type::FKL) continue;
ev->central().weight *= invGeV2_to_pb;
ev->central().weight /= config.events;
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 << std::endl;
assert(std::abs(xs - xs_ref) < 3*xs_err);
assert(xs_err < 0.05*xs);
return EXIT_SUCCESS;
}
diff --git a/FixedOrderGen/t/h_5j.cc b/FixedOrderGen/t/h_5j.cc
index 851359a..0e0948f 100644
--- a/FixedOrderGen/t/h_5j.cc
+++ b/FixedOrderGen/t/h_5j.cc
@@ -1,72 +1,72 @@
/**
* This is a regression test
* the reference cross section has not been checked against any other program
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#ifdef NDEBUG
#undef NDEBUG
#endif
#include <algorithm>
#include <cassert>
#include <cmath>
#include <iostream>
#include "config.hh"
#include "EventGenerator.hh"
#include "HEJ/Ranlux64.hh"
#include "HEJ/Event.hh"
#include "HEJ/MatrixElement.hh"
#include "HEJ/PDF.hh"
using namespace HEJFOG;
int main(){
constexpr double invGeV2_to_pb = 389379292.;
constexpr double xs_ref = 0.252273; // +- 0.00657742
//calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20
auto config = load_config("config_h_2j.yml");
config.process.njets = 5;
- HEJ::Ranlux64 ran{};
+ std::shared_ptr<HEJ::RNG> ran{std::make_shared<HEJ::Ranlux64>()};
HEJFOG::EventGenerator generator{
config.process,
config.beam,
HEJ::ScaleGenerator{
config.scales.base,
config.scales.factors,
config.scales.max_ratio
},
config.jets,
config.pdf_id,
config.subleading_fraction,
config.subleading_channels,
config.particle_decays,
config.Higgs_coupling,
config.ew_parameters,
ran
};
double xs = 0., xs_err = 0.;
for (size_t trials = 0; trials < config.events; ++trials){
auto ev = generator.gen_event();
if(generator.status() != good) continue;
assert(ev);
ev->central().weight *= invGeV2_to_pb;
ev->central().weight /= config.events;
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 << std::endl;
assert(std::abs(xs - xs_ref) < 3*xs_err);
assert(xs_err < 0.06*xs);
return EXIT_SUCCESS;
}
diff --git a/include/HEJ/EventReweighter.hh b/include/HEJ/EventReweighter.hh
index ed5a953..87a0279 100644
--- a/include/HEJ/EventReweighter.hh
+++ b/include/HEJ/EventReweighter.hh
@@ -1,201 +1,197 @@
/** \file
* \brief Declares the EventReweighter class
*
* EventReweighter is the main class used within HEJ 2. It reweights the
* resummation events.
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <array>
-#include <functional>
-#include <utility>
+#include <memory>
#include <vector>
#include "HEJ/Config.hh"
#include "HEJ/event_types.hh"
#include "HEJ/MatrixElement.hh"
#include "HEJ/Parameters.hh"
#include "HEJ/PDF.hh"
#include "HEJ/PDG_codes.hh"
-#include "HEJ/RNG.hh"
#include "HEJ/ScaleFunction.hh"
#include "HEJ/StatusCode.hh"
namespace LHEF {
class HEPRUP;
}
namespace HEJ {
class Event;
+ class RNG;
//! Beam parameters
/**
* Currently, only symmetric beams are supported,
* so there is a single beam energy.
*/
struct Beam{
double E; /**< Beam energy */
std::array<ParticleID, 2> type; /**< Beam particles */
};
//! Main class for reweighting events in HEJ.
class EventReweighter{
using EventType = event_type::EventType;
public:
EventReweighter(
Beam beam, /**< Beam Energy */
int pdf_id, /**< PDF ID */
ScaleGenerator scale_gen, /**< Scale settings */
EventReweighterConfig conf, /**< Configuration parameters */
- HEJ::RNG & ran /**< Random number generator */
+ std::shared_ptr<RNG> ran /**< Random number generator */
);
EventReweighter(
LHEF::HEPRUP const & heprup, /**< LHEF event header */
ScaleGenerator scale_gen, /**< Scale settings */
EventReweighterConfig conf, /**< Configuration parameters */
- HEJ::RNG & ran /**< Random number generator */
+ std::shared_ptr<RNG> ran /**< Random number generator */
);
//! Get the used pdf
PDF const & pdf() const;
//! Get event treatment
EventTreatment treatment(EventType type) const;
//! Generate resummation events for a given fixed-order event
/**
* @param ev Fixed-order event corresponding
* to the resummation events
* @param num_events Number of trial resummation configurations.
* @returns A vector of resummation events.
*
* The result vector depends on the type of the input event and the
* treatment of different types as specified in the constructor:
*
* \ref reweight The result vector contains between
* 0 and num_events resummation events.
*
* \ref keep If the input event passes the resummation jet cuts
* the result vector contains one event. Otherwise it is empty.
*
* \ref discard The result vector is empty
*/
std::vector<Event> reweight(
Event const & ev,
int num_events
);
//! Gives all StatusCodes of the last reweight()
/**
* Each StatusCode corresponds to one tried generation. Only good
* StatusCodes generated an event.
*/
std::vector<StatusCode> const & status() const {
return status_;
}
private:
template<typename... T>
PDF const & pdf(T&& ...);
/** \internal
* \brief main generation/reweighting function:
* generate phase space points and divide out Born factors
*/
std::vector<Event> gen_res_events(
Event const & ev, size_t num_events
);
std::vector<Event> rescale(
Event const & Born_ev, std::vector<Event> events
) const;
/** \internal
* \brief Do the Jets pass the resummation Cuts?
*
* @param ev Event in Question
* @returns 0 or 1 depending on if ev passes Jet Cuts
*/
bool jets_pass_resummation_cuts(Event const & ev) const;
/** \internal
* \brief pdf_factors Function
*
* @param ev Event in Question
* @returns EventFactor due to PDFs
*
* Calculates the Central value and the variation due
* to the PDF choice made.
*/
Weights pdf_factors(Event const & ev) const;
/** \internal
* \brief matrix_elements Function
*
* @param ev Event in question
* @returns EventFactor due to MatrixElements
*
* Calculates the Central value and the variation due
* to the Matrix Element.
*/
Weights matrix_elements(Event const & ev) const;
/** \internal
* \brief Scale-dependent part of fixed-order matrix element
*
* @param ev Event in question
* @returns EventFactor scale variation due to FO-ME.
*
* This is only called to compute the scale variation for events where
* we don't do resummation (e.g. non-FKL).
* Since at tree level the scale dependence is just due to alpha_s,
* it is enough to return the alpha_s(mur) factors in the matrix element.
* The rest drops out in the ratio of (output event ME)/(input event ME),
* so we never have to compute it.
*/
Weights fixed_order_scale_ME(Event const & ev) const;
/** \internal
* \brief Computes the tree level matrix element
*
* @param ev Event in Question
* @returns HEJ approximation to Tree level Matrix Element
*
* This computes the HEJ approximation to the tree level FO
* Matrix element which is used within the LO weighting process.
*/
double tree_matrix_element(Event const & ev) const;
//! \internal General parameters
EventReweighterConfig param_;
//! \internal Beam energy
double E_beam_;
//! \internal PDF
PDF pdf_;
//! \internal Object to calculate the square of the matrix element
MatrixElement MEt2_;
//! \internal Object to calculate event renormalisation and factorisation scales
ScaleGenerator scale_gen_;
- /** \internal random number generator
- *
- * \note We use a reference_wrapper so that EventReweighter objects can
- * still be copied (which would be impossible with a reference).
- */
- std::reference_wrapper<HEJ::RNG> ran_;
+ //! \internal random number generator
+ std::shared_ptr<RNG> ran_;
+ //! \internal StatusCode of each attempt
std::vector<StatusCode> status_;
};
template<typename... T>
PDF const & EventReweighter::pdf(T&&... t){
return pdf_ = PDF{std::forward<T>(t)...};
}
} // namespace HEJ
diff --git a/include/HEJ/JetSplitter.hh b/include/HEJ/JetSplitter.hh
index f634d7c..bf0fabd 100644
--- a/include/HEJ/JetSplitter.hh
+++ b/include/HEJ/JetSplitter.hh
@@ -1,77 +1,70 @@
/**
* \file
* \brief Declaration of the JetSplitter class
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
-#include <functional>
+#include <memory>
#include <vector>
#include "fastjet/JetDefinition.hh"
-#include "HEJ/RNG.hh"
-
namespace fastjet {
class PseudoJet;
}
namespace HEJ {
+ class RNG;
+
//! Class to split jets into their constituents
class JetSplitter {
public:
struct SplitResult {
std::vector<fastjet::PseudoJet> constituents;
double weight;
};
//! Constructor
/**
- * @param jet_def Jet definition
- * @param min_pt Minimum jet transverse momentum
- * @param ran Random number generator
+ * @param jet_def Jet definition
+ * @param min_pt Minimum jet transverse momentum
*/
- JetSplitter(
- fastjet::JetDefinition jet_def, double min_pt,
- HEJ::RNG & ran
- ):
- R_{jet_def.R()},
- min_jet_pt_{min_pt},
- jet_def_{jet_def},
- ran_{ran}
- {}
+ JetSplitter(fastjet::JetDefinition jet_def, double min_pt);
//! Split a get into constituents
/**
- * @param j2split Jet to be split
- * @param ncons Number of constituents
- * @returns The constituent momenta
- * together with the associated weight
+ * @param j2split Jet to be split
+ * @param ncons Number of constituents
+ * @param ran Random number generator
+ * @returns The constituent momenta together with the associated
+ * weight
*/
- SplitResult split(fastjet::PseudoJet const & j2split, int ncons) const;
+ SplitResult split(
+ fastjet::PseudoJet const & j2split, int ncons, RNG & ran
+ ) const;
//! Maximum distance of constituents to jet axis
static constexpr double R_factor = 5./3.;
private:
//! \internal split jet into two partons
- SplitResult Split2(fastjet::PseudoJet const & j2split) const;
+ SplitResult Split2(fastjet::PseudoJet const & j2split, RNG & ran) const;
/** \internal
* @brief sample y-phi distance to jet pt axis for a jet splitting into two
* partons
*
* @param wt Multiplied by the weight of the sampling point
+ * @param ran Random number generator
* @returns The distance in units of the jet radius
*/
- double sample_distance_2p(double & wt) const;
+ double sample_distance_2p(double & wt, RNG & ran) const;
- double R_;
double min_jet_pt_;
fastjet::JetDefinition jet_def_;
- std::reference_wrapper<HEJ::RNG> ran_;
};
} // namespace HEJ
diff --git a/include/HEJ/PhaseSpacePoint.hh b/include/HEJ/PhaseSpacePoint.hh
index 19ad82e..389f803 100644
--- a/include/HEJ/PhaseSpacePoint.hh
+++ b/include/HEJ/PhaseSpacePoint.hh
@@ -1,197 +1,198 @@
/** \file
* \brief Contains the PhaseSpacePoint Class
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <array>
-#include <functional>
#include <unordered_map>
#include <vector>
#include "HEJ/Config.hh"
#include "HEJ/Particle.hh"
-#include "HEJ/RNG.hh"
#include "HEJ/StatusCode.hh"
namespace HEJ{
class Event;
+ class RNG;
//! A point in resummation phase space
class PhaseSpacePoint{
public:
//! No default PhaseSpacePoint Constructor
PhaseSpacePoint() = delete;
//! PhaseSpacePoint Constructor
/**
* @param ev Clustered Jet Event
* @param conf Configuration parameters
* @param ran Random number generator
*/
PhaseSpacePoint(
Event const & ev,
PhaseSpacePointConfig conf,
RNG & ran
);
//! Get phase space point weight
double weight() const{
return weight_;
}
//! Access incoming particles
std::array<Particle, 2> const & incoming() const{
return incoming_;
}
//! Access outgoing particles
std::vector<Particle> const & outgoing() const{
return outgoing_;
}
//! Particle decays
/**
* The key in the returned map corresponds to the index in the
* vector returned by outgoing()
*/
std::unordered_map<size_t, std::vector<Particle>> const & decays() const{
return decays_;
}
//! Status code of generation
StatusCode status() const{
return status_;
}
static constexpr int ng_max = 1000; //!< maximum number of extra gluons
private:
//! /internal returns the clustered jets sorted in rapidity
std::vector<fastjet::PseudoJet> cluster_jets(
std::vector<fastjet::PseudoJet> const & partons
) const;
bool pass_resummation_cuts(
std::vector<fastjet::PseudoJet> const & jets
) const;
bool pass_extremal_cuts(
fastjet::PseudoJet const & ext_parton,
fastjet::PseudoJet const & jet
) const;
- int sample_ng(std::vector<fastjet::PseudoJet> const & Born_jets);
- int sample_ng_jets(int ng, std::vector<fastjet::PseudoJet> const & Born_jets);
+ int sample_ng(std::vector<fastjet::PseudoJet> const & Born_jets, RNG & ran);
+ int sample_ng_jets(
+ int ng, std::vector<fastjet::PseudoJet> const & Born_jets, RNG & ran
+ );
double probability_in_jet(
std::vector<fastjet::PseudoJet> const & Born_jets
) const;
std::vector<fastjet::PseudoJet> gen_non_jet(
int ng_non_jet,
- double ptmin, double ptmax
+ double ptmin, double ptmax,
+ RNG & ran
);
void rescale_qqx_rapidities(
std::vector<fastjet::PseudoJet> & out_partons,
std::vector<fastjet::PseudoJet> const & jets,
const double ymin1, const double ymax2,
const int qqxbackjet
);
void rescale_rapidities(
std::vector<fastjet::PseudoJet> & partons,
double ymin, double ymax
);
std::vector<fastjet::PseudoJet> reshuffle(
std::vector<fastjet::PseudoJet> const & Born_jets,
fastjet::PseudoJet const & q
);
/** \interal final jet test
* - number of jets must match Born kinematics
* - no partons designated as nonjet may end up inside jets
* - all other outgoing partons *must* end up inside jets
* - the extremal (in rapidity) partons must be inside the extremal jets
* - rapidities must be the same (by construction)
*/
bool jets_ok(
std::vector<fastjet::PseudoJet> const & Born_jets,
std::vector<fastjet::PseudoJet> const & partons
) const;
void reconstruct_incoming(std::array<Particle, 2> const & Born_incoming);
double phase_space_normalisation(
int num_Born_jets,
int num_res_partons
) const;
/** \interal Distribute gluon in jet
* @param jets jets to distribute gluon in
* @param ng_jets number of gluons
* @param qqxbackjet position of first (backwards) qqx jet
*
* relies on JetSplitter
*/
std::vector<fastjet::PseudoJet> split(
std::vector<fastjet::PseudoJet> const & jets,
- int ng_jets, size_t qqxbackjet
+ int ng_jets, size_t qqxbackjet, RNG & ran
);
std::vector<int> distribute_jet_partons(
- int ng_jets, std::vector<fastjet::PseudoJet> const & jets
+ int ng_jets, std::vector<fastjet::PseudoJet> const & jets, RNG & ran
);
std::vector<fastjet::PseudoJet> split(
std::vector<fastjet::PseudoJet> const & jets,
std::vector<int> const & np_in_jet,
- size_t qqxbackjet
+ size_t qqxbackjet,
+ RNG & ran
);
bool split_preserved_jets(
std::vector<fastjet::PseudoJet> const & jets,
std::vector<fastjet::PseudoJet> const & jet_partons
) const;
template<class Particle>
Particle const & most_backward_FKL(
std::vector<Particle> const & partons
) const;
template<class Particle>
Particle const & most_forward_FKL(
std::vector<Particle> const & partons
) const;
template<class Particle>
Particle & most_backward_FKL(std::vector<Particle> & partons) const;
template<class Particle>
Particle & most_forward_FKL(std::vector<Particle> & partons) const;
bool extremal_ok(
std::vector<fastjet::PseudoJet> const & partons
) const;
/** \internal
* Assigns PDG IDs to outgoing partons, i.e. labels them as quarks
*
* \note This function assumes outgoing_ to be pure partonic when called,
* i.e. A/W/Z/h bosons should _not be set_ at this stage
*/
void label_quarks(Event const & event);
/** \internal
* This function will label the qqx pair in a qqx event back to their
* original types from the input event.
*/
void label_qqx(Event const & event);
void copy_AWZH_boson_from(Event const & event);
bool momentum_conserved() const;
bool contains_idx(
fastjet::PseudoJet const & jet, fastjet::PseudoJet const & parton
) const;
bool unob_, unof_, qqxb_, qqxf_, qqxmid_;
double weight_;
PhaseSpacePointConfig param_;
std::array<Particle, 2> incoming_;
std::vector<Particle> outgoing_;
//! \internal Particle decays in the format {outgoing index, decay products}
std::unordered_map<size_t, std::vector<Particle>> decays_;
- std::reference_wrapper<HEJ::RNG> ran_;
-
StatusCode status_;
};
} // namespace HEJ
diff --git a/include/HEJ/make_RNG.hh b/include/HEJ/make_RNG.hh
index ec7c026..814b6c2 100644
--- a/include/HEJ/make_RNG.hh
+++ b/include/HEJ/make_RNG.hh
@@ -1,32 +1,32 @@
/** \file
* \brief Declares a factory function for random number generators
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <memory>
#include <string>
#include "HEJ/optional.hh"
#include "HEJ/RNG.hh"
namespace HEJ {
//! Factory function for random number generators
/**
* @param name Name of the random number generator
* @param seed Optional seed
* @returns A pointer to an instance of a random number generator
*
* At present, name should be one of "ranlux64" or "mixmax" (case insensitive).
* The interpretation of the seed depends on the random number generator.
* For ranlux64, it is the name of a seed file. For mixmax it should be a
* string convertible to a long integer.
*/
- std::unique_ptr<HEJ::RNG> make_RNG(
+ std::unique_ptr<RNG> make_RNG(
std::string const & name,
optional<std::string> const & seed
);
}
diff --git a/src/EventReweighter.cc b/src/EventReweighter.cc
index 08fddb0..bf28520 100644
--- a/src/EventReweighter.cc
+++ b/src/EventReweighter.cc
@@ -1,282 +1,284 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/EventReweighter.hh"
#include <algorithm>
#include <assert.h>
#include <limits>
#include <math.h>
#include <stddef.h>
#include <string>
#include <unordered_map>
+#include <utility>
#include "fastjet/ClusterSequence.hh"
#include "LHEF/LHEF.h"
#include "HEJ/Event.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/Particle.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/PhaseSpacePoint.hh"
+#include "HEJ/RNG.hh"
namespace HEJ{
- using EventType = event_type::EventType;
-
namespace {
static_assert(
std::numeric_limits<double>::has_quiet_NaN,
"no quiet NaN for double"
);
constexpr double NaN = std::numeric_limits<double>::quiet_NaN();
Event::EventData to_EventData(PhaseSpacePoint const & psp){
Event::EventData result;
result.incoming=psp.incoming();
assert(result.incoming.size() == 2);
result.outgoing=psp.outgoing();
// technically Event::EventData doesn't have to be sorted,
// but PhaseSpacePoint should be anyway
assert(
std::is_sorted(
begin(result.outgoing), end(result.outgoing),
rapidity_less{}
)
);
assert(result.outgoing.size() >= 2);
result.decays = psp.decays();
result.parameters.central = {NaN, NaN, psp.weight()};
return result;
}
} // namespace anonymous
EventReweighter::EventReweighter(
LHEF::HEPRUP const & heprup,
ScaleGenerator scale_gen,
EventReweighterConfig conf,
- HEJ::RNG & ran
+ std::shared_ptr<RNG> ran
):
EventReweighter{
HEJ::Beam{
heprup.EBMUP.first,
{{
static_cast<HEJ::ParticleID>(heprup.IDBMUP.first),
static_cast<HEJ::ParticleID>(heprup.IDBMUP.second)
}}
},
heprup.PDFSUP.first,
std::move(scale_gen),
std::move(conf),
- ran
+ std::move(ran)
}
{
if(heprup.EBMUP.second != E_beam_){
throw std::invalid_argument(
"asymmetric beam: " + std::to_string(E_beam_)
+ " ---> <--- " + std::to_string(heprup.EBMUP.second)
);
};
if(heprup.PDFSUP.second != pdf_.id()){
throw std::invalid_argument(
"conflicting PDF ids: " + std::to_string(pdf_.id())
+ " vs. " + std::to_string(heprup.PDFSUP.second)
);
}
}
EventReweighter::EventReweighter(
Beam beam,
int pdf_id,
ScaleGenerator scale_gen,
EventReweighterConfig conf,
- HEJ::RNG & ran
+ std::shared_ptr<RNG> ran
):
param_{std::move(conf)},
E_beam_{beam.E},
pdf_{pdf_id, beam.type.front(), beam.type.back()},
MEt2_{
[this](double mu){ return pdf_.Halphas(mu); },
param_.ME_config
},
- scale_gen_(std::move(scale_gen)),
- ran_{ran}
- {}
+ scale_gen_{std::move(scale_gen)},
+ ran_{std::move(ran)}
+ {
+ assert(ran_);
+ }
PDF const & EventReweighter::pdf() const{
return pdf_;
}
std::vector<Event> EventReweighter::reweight(
Event const & input_ev, int num_events
){
auto res_events{ gen_res_events(input_ev, num_events) };
if(res_events.empty()) return {};
for(auto & event: res_events) event = scale_gen_(std::move(event));
return rescale(input_ev, std::move(res_events));
}
EventTreatment EventReweighter::treatment(EventType type) const {
return param_.treat.at(type);
}
std::vector<Event> EventReweighter::gen_res_events(
Event const & ev,
size_t phase_space_points
){
assert(ev.variations().empty());
status_.clear();
switch(treatment(ev.type())){
case EventTreatment::discard: {
status_.emplace_back(StatusCode::discard);
return {};
}
case EventTreatment::keep:
if(! jets_pass_resummation_cuts(ev)) {
status_.emplace_back(StatusCode::failed_resummation_cuts);
return {};
}
else {
status_.emplace_back(StatusCode::good);
return {ev};
}
default:;
}
const double Born_shat = shat(ev);
std::vector<Event> resummation_events;
status_.reserve(phase_space_points);
for(size_t psp_number = 0; psp_number < phase_space_points; ++psp_number){
- PhaseSpacePoint psp{ev, param_.psp_config, ran_};
+ PhaseSpacePoint psp{ev, param_.psp_config, *ran_};
status_.emplace_back(psp.status());
assert(psp.status() != StatusCode::unspecified);
if(psp.status() != StatusCode::good) continue;
assert(psp.weight() != 0.);
if(psp.incoming()[0].E() > E_beam_ || psp.incoming()[1].E() > E_beam_) {
status_.back() = StatusCode::too_much_energy;
continue;
}
resummation_events.emplace_back(
to_EventData( std::move(psp) ).cluster(
param_.jet_param().def, param_.jet_param().min_pt
)
);
auto & new_event = resummation_events.back();
assert( new_event.valid_hej_state(
param_.psp_config.max_ext_soft_pt_fraction,
param_.psp_config.min_extparton_pt ) );
if( new_event.type() != ev.type() )
throw std::logic_error{"Resummation Event does not match Born event"};
- new_event.generate_colours(ran_);
+ new_event.generate_colours(*ran_);
assert(new_event.variations().empty());
new_event.central().mur = ev.central().mur;
new_event.central().muf = ev.central().muf;
const double resum_shat = shat(new_event);
new_event.central().weight *= ev.central().weight*Born_shat*Born_shat/
(phase_space_points*resum_shat*resum_shat);
}
return resummation_events;
}
std::vector<Event> EventReweighter::rescale(
Event const & Born_ev,
std::vector<Event> events
) const{
const double Born_pdf = pdf_factors(Born_ev).central;
const double Born_ME = tree_matrix_element(Born_ev);
for(auto & cur_event: events){
const auto pdf = pdf_factors(cur_event);
assert(pdf.variations.size() == cur_event.variations().size());
const auto ME = matrix_elements(cur_event);
assert(ME.variations.size() == cur_event.variations().size());
cur_event.parameters() *= pdf*ME/(Born_pdf*Born_ME);
}
return events;
}
bool EventReweighter::jets_pass_resummation_cuts(
Event const & ev
) const{
const auto out_as_PseudoJet = to_PseudoJet(filter_partons(ev.outgoing()));
fastjet::ClusterSequence cs{out_as_PseudoJet, param_.jet_param().def};
return cs.inclusive_jets(param_.jet_param().min_pt).size() == ev.jets().size();
}
Weights EventReweighter::pdf_factors(Event const & ev) const{
auto const & a = ev.incoming().front();
auto const & b = ev.incoming().back();
const double xa = a.p.e()/E_beam_;
const double xb = b.p.e()/E_beam_;
Weights result;
std::unordered_map<double, double> known_pdf;
result.central =
pdf_.pdfpt(0,xa,ev.central().muf,a.type)*
pdf_.pdfpt(1,xb,ev.central().muf,b.type);
known_pdf.emplace(ev.central().muf, result.central);
result.variations.reserve(ev.variations().size());
for(auto const & ev_param: ev.variations()){
const double muf = ev_param.muf;
auto cur_pdf = known_pdf.find(muf);
if(cur_pdf == known_pdf.end()){
cur_pdf = known_pdf.emplace(
muf,
pdf_.pdfpt(0,xa,muf,a.type)*pdf_.pdfpt(1,xb,muf,b.type)
).first;
}
result.variations.emplace_back(cur_pdf->second);
}
assert(result.variations.size() == ev.variations().size());
return result;
}
Weights
EventReweighter::matrix_elements(Event const & ev) const{
assert(param_.treat.count(ev.type()) > 0);
if(param_.treat.find(ev.type())->second == EventTreatment::keep){
return fixed_order_scale_ME(ev);
}
return MEt2_(ev);
}
double EventReweighter::tree_matrix_element(Event const & ev) const{
assert(ev.variations().empty());
assert(param_.treat.count(ev.type()) > 0);
if(param_.treat.find(ev.type())->second == EventTreatment::keep){
return fixed_order_scale_ME(ev).central;
}
return MEt2_.tree(ev).central;
}
Weights
EventReweighter::fixed_order_scale_ME(Event const & ev) const{
int alpha_s_power = 0;
for(auto const & part: ev.outgoing()){
if(is_parton(part))
++alpha_s_power;
else if(part.type == pid::Higgs) {
alpha_s_power += 2;
}
// nothing to do for other uncoloured particles
}
Weights result;
result.central = pow(pdf_.Halphas(ev.central().mur), alpha_s_power);
for(auto const & var: ev.variations()){
result.variations.emplace_back(
pow(pdf_.Halphas(var.mur), alpha_s_power)
);
}
return result;
}
} // namespace HEJ
diff --git a/src/JetSplitter.cc b/src/JetSplitter.cc
index 60003ce..96a9162 100644
--- a/src/JetSplitter.cc
+++ b/src/JetSplitter.cc
@@ -1,178 +1,188 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/JetSplitter.hh"
#include <array>
#include <assert.h>
#include <numeric>
+#include <utility>
#include "fastjet/ClusterSequence.hh"
#include "fastjet/PseudoJet.hh"
#include "HEJ/Constants.hh"
#include "HEJ/exceptions.hh"
+#include "HEJ/RNG.hh"
namespace HEJ {
namespace{
constexpr double ccut=HEJ::CMINPT; // min parton pt
template<class Iterator>
bool same_pt_and_rapidity(
Iterator begin, Iterator end,
fastjet::PseudoJet const & jet
){
constexpr double ep = 1e-2;
const fastjet::PseudoJet reconstructed_jet = std::accumulate(
begin, end, fastjet::PseudoJet{}
);
return
(std::abs(reconstructed_jet.pt() - jet.pt()) < ep)
&& (std::abs(reconstructed_jet.rapidity() - jet.rapidity()) < ep)
;
}
bool all_in_one_jet(
std::vector<fastjet::PseudoJet> const & partons,
fastjet::JetDefinition jet_def, double min_jet_pt
){
fastjet::ClusterSequence ev(partons, jet_def);
const std::vector<fastjet::PseudoJet> testjet = ev.inclusive_jets(min_jet_pt);
return testjet.size() == 1u
&& testjet[0].constituents().size() == partons.size();
}
}
+ JetSplitter::JetSplitter(
+ fastjet::JetDefinition jet_def, double min_pt
+ ):
+ min_jet_pt_{min_pt},
+ jet_def_{std::move(jet_def)}
+ {}
using SplitResult = JetSplitter::SplitResult;
SplitResult JetSplitter::split(
- fastjet::PseudoJet const & j2split, int ncons
+ fastjet::PseudoJet const & j2split, int ncons, RNG & ran
) const{
if(ncons <= 0) {
throw std::invalid_argument{
"number of requested jet constituents less than 1"
};
}
double swt = 1.;
std::vector<fastjet::PseudoJet> jcons;
if(ncons == 1){
jcons.emplace_back(j2split);
jcons.back().set_user_index(0);
return {jcons, swt};
}
if(ncons == 2){
- return Split2(j2split);
+ return Split2(j2split, ran);
}
- const double R_max = R_factor*R_;
+ const double R_max = R_factor*jet_def_.R();
assert(R_max < M_PI);
double pt_remaining = j2split.pt();
const double phi_jet = j2split.phi();
const double y_jet = j2split.rapidity();
for(int i = 0; i < ncons - 1; ++i){
/**
* Generate rapidity and azimuthal angle with a distance
* R = sqrt(delta_y^2 + delta_phi^2) < R_max
* from the jet centre
*/
- const double R = R_max*ran_.get().flat();
- const double theta = 2*M_PI*ran_.get().flat();
+ const double R = R_max*ran.flat();
+ const double theta = 2*M_PI*ran.flat();
const double delta_phi = R*cos(theta);
const double delta_y = R*sin(theta);
/**
* Generate pt such that the total contribution of all partons
* along the jet pt axis does not exceed the jet pt
*/
const double pt_max = pt_remaining/cos(delta_phi);
assert(pt_max > 0);
if(pt_max < ccut) return {}; // no pt remaining for this parton
- const double pt = (pt_max - ccut)*ran_.get().flat() + ccut;
+ const double pt = (pt_max - ccut)*ran.flat() + ccut;
pt_remaining -= pt*cos(delta_phi);
jcons.emplace_back(
pt*cos(phi_jet + delta_phi), pt*sin(phi_jet + delta_phi),
pt*sinh(y_jet + delta_y), pt*cosh(y_jet + delta_y)
);
jcons.back().set_user_index(i);
swt *= 2*M_PI*R*R_max*pt*(pt_max - ccut);
}
const fastjet::PseudoJet p_total = std::accumulate(
jcons.begin(), jcons.end(), fastjet::PseudoJet{}
);
// Calculate the pt of the last parton
const double last_px = j2split.px() - p_total.px();
const double last_py = j2split.py() - p_total.py();
const double last_pt = sqrt(last_px*last_px + last_py*last_py);
if(last_pt < ccut) return {};
// Calculate the rapidity of the last parton using the requirement that the
// new jet must have the same rapidity as the LO jet.
const double exp_2y_jet = (j2split.e() + j2split.pz())/(j2split.e() - j2split.pz());
const double bb = (p_total.e()+p_total.pz()) - exp_2y_jet*(p_total.e()-p_total.pz());
const double lasty = log((-bb+sqrt(bb*bb+4.*exp_2y_jet*last_pt*last_pt))/(2.*last_pt));
jcons.emplace_back(
last_px, last_py, last_pt*sinh(lasty), last_pt*cosh(lasty)
);
jcons.back().set_user_index(ncons-1);
assert(same_pt_and_rapidity(begin(jcons), end(jcons), j2split));
// Test that the last parton is not too far away from the jet centre.
if (jcons.back().delta_R(j2split) > R_max) return {};
if(! all_in_one_jet(jcons, jet_def_, min_jet_pt_)) return {};
return {jcons, swt};
}
- double JetSplitter::sample_distance_2p(double & wt) const{
+ double JetSplitter::sample_distance_2p(double & wt, RNG & ran) const{
static constexpr double x_small = 0.1;
static constexpr double p_small = 0.4;
- const double pR = ran_.get().flat();
+ const double pR = ran.flat();
if(pR < p_small){
wt *= x_small/p_small;
return x_small/p_small*pR;
}
wt *= (1-x_small)/(1-p_small);
return (1-x_small)/(1-p_small)*(pR-p_small) + x_small;
}
- SplitResult JetSplitter::Split2(fastjet::PseudoJet const & j2split) const{
+ SplitResult JetSplitter::Split2(
+ fastjet::PseudoJet const & j2split, RNG & ran
+ ) const{
static constexpr size_t ncons = 2;
std::vector<fastjet::PseudoJet> jcons(ncons);
std::array<double, ncons> R, phi, y, pt;
double wt = 1;
- const double theta = 2*M_PI*ran_.get().flat(); // angle in y-phi plane
+ const double theta = 2*M_PI*ran.flat(); // angle in y-phi plane
// empiric observation: we are always within the jet radius
- R[0] = sample_distance_2p(wt)*R_;
- R[1] = -sample_distance_2p(wt)*R_;
+ R[0] = sample_distance_2p(wt, ran)*jet_def_.R();
+ R[1] = -sample_distance_2p(wt, ran)*jet_def_.R();
for(size_t i = 0; i <= 1; ++i){
phi[i] = j2split.phi() + R[i]*cos(theta);
y[i] = j2split.rapidity() + R[i]*sin(theta);
}
for(size_t i = 0; i <= 1; ++i){
pt[i] = (j2split.py() - tan(phi[1-i])*j2split.px())/
(sin(phi[i]) - tan(phi[1-i])*cos(phi[i]));
if(pt[i] < ccut) return {};
jcons[i].reset_PtYPhiM(pt[i], y[i], phi[i]);
jcons[i].set_user_index(i);
}
assert(same_pt_and_rapidity(begin(jcons), end(jcons), j2split));
if(! all_in_one_jet(jcons, jet_def_, min_jet_pt_)) return {};
- wt *= 2*M_PI*pt[0]*R[0]*R_*R_;
+ wt *= 2*M_PI*pt[0]*R[0]*jet_def_.R()*jet_def_.R();
// from transformation of delta(R[1] - ...) to delta(pt[0] - ...)
const double dphi0 = phi[0] - j2split.phi();
const double ptJ = j2split.pt();
const double jacobian = cos(theta)*pt[1]*pt[1]/(ptJ*sin(dphi0));
return {jcons, jacobian*wt};
}
}
diff --git a/src/PhaseSpacePoint.cc b/src/PhaseSpacePoint.cc
index 726c0b7..9c5b669 100644
--- a/src/PhaseSpacePoint.cc
+++ b/src/PhaseSpacePoint.cc
@@ -1,810 +1,814 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/PhaseSpacePoint.hh"
#include <algorithm>
#include <assert.h>
#include <numeric>
#include <random>
#include "fastjet/ClusterSequence.hh"
#include "HEJ/Constants.hh"
#include "HEJ/Event.hh"
#include "HEJ/JetSplitter.hh"
#include "HEJ/kinematics.hh"
#include "HEJ/resummation_jet.hh"
#include "HEJ/utility.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/event_types.hh"
namespace HEJ{
namespace {
constexpr int max_jet_user_idx = PhaseSpacePoint::ng_max;
bool is_nonjet_parton(fastjet::PseudoJet const & parton){
assert(parton.user_index() != -1);
return parton.user_index() > max_jet_user_idx;
}
bool is_jet_parton(fastjet::PseudoJet const & parton){
assert(parton.user_index() != -1);
return parton.user_index() <= max_jet_user_idx;
}
// user indices for partons with extremal rapidity
constexpr int qqxmid1_idx = -9;
constexpr int qqxmid2_idx = -8;
constexpr int qqxb_idx = -7;
constexpr int qqxf_idx = -6;
constexpr int unob_idx = -5;
constexpr int unof_idx = -4;
constexpr int backward_FKL_idx = -3;
constexpr int forward_FKL_idx = -2;
}
namespace {
double estimate_ng_mean(std::vector<fastjet::PseudoJet> const & Born_jets){
const double delta_y =
Born_jets.back().rapidity() - Born_jets.front().rapidity();
assert(delta_y > 0);
// Formula derived from fit in arXiv:1805.04446 (see Fig. 2)
return 0.975052*delta_y;
}
}
std::vector<fastjet::PseudoJet> PhaseSpacePoint::cluster_jets(
std::vector<fastjet::PseudoJet> const & partons
) const{
fastjet::ClusterSequence cs(partons, param_.jet_param.def);
return sorted_by_rapidity(cs.inclusive_jets(param_.jet_param.min_pt));
}
bool PhaseSpacePoint::pass_resummation_cuts(
std::vector<fastjet::PseudoJet> const & jets
) const{
return cluster_jets(jets).size() == jets.size();
}
- int PhaseSpacePoint::sample_ng(std::vector<fastjet::PseudoJet> const & Born_jets){
+ int PhaseSpacePoint::sample_ng(
+ std::vector<fastjet::PseudoJet> const & Born_jets, RNG & ran
+ ){
const double ng_mean = estimate_ng_mean(Born_jets);
std::poisson_distribution<int> dist(ng_mean);
- const int ng = dist(ran_.get());
+ const int ng = dist(ran);
assert(ng >= 0);
assert(ng < ng_max);
weight_ *= std::tgamma(ng + 1)*std::exp(ng_mean)*std::pow(ng_mean, -ng);
return ng;
}
void PhaseSpacePoint::copy_AWZH_boson_from(Event const & event){
auto const & from = event.outgoing();
const auto AWZH_boson = std::find_if(
begin(from), end(from),
[](Particle const & p){ return is_AWZH_boson(p); }
);
if(AWZH_boson == end(from)) return;
auto insertion_point = std::lower_bound(
begin(outgoing_), end(outgoing_), *AWZH_boson, rapidity_less{}
);
outgoing_.insert(insertion_point, *AWZH_boson);
// copy decay products
const int idx = std::distance(begin(from), AWZH_boson);
assert(idx >= 0);
const auto decay_it = event.decays().find(idx);
if(decay_it != end(event.decays())){
const int new_idx = std::distance(begin(outgoing_), insertion_point);
assert(new_idx >= 0);
assert(outgoing_[new_idx].type == AWZH_boson->type);
decays_.emplace(new_idx, decay_it->second);
}
assert(std::is_sorted(begin(outgoing_), end(outgoing_), rapidity_less{}));
}
namespace {
auto get_first_anyquark_emission(Event const & ev) {
// find born quarks (ignore extremal partons)
auto const firstquark = std::find_if(
std::next(ev.begin_partons()), std::prev(ev.end_partons(), 2),
[](Particle const & s){ return (is_anyquark(s)); }
);
// assert that it is a q-q_bar pair.
assert(std::distance(firstquark, ev.end_partons()) != 2);
assert(
( is_quark(*firstquark) && is_antiquark(*std::next(firstquark)) )
|| ( is_antiquark(*firstquark) && is_quark(*std::next(firstquark)) )
);
return firstquark;
}
//! returns index of most backward q-qbar jet
template<class Iterator>
int get_back_quark_jet(Event const & ev, Iterator firstquark){
// find jets at FO corresponding to the quarks
// technically this isn't necessary for LO
std::vector<int> const born_indices{ ev.particle_jet_indices() };
const auto firstquark_idx = std::distance(ev.begin_partons(), firstquark);
int const firstjet_idx = born_indices[firstquark_idx];
assert(firstjet_idx>0);
assert( born_indices[firstquark_idx+1] == firstjet_idx+1 );
return firstjet_idx;
}
//! returns index of most backward q-qbar jet
int getBackQuarkJet(Event const & ev){
const auto firstquark = get_first_anyquark_emission(ev);
return get_back_quark_jet(ev, firstquark);
}
template<class ConstIterator, class Iterator>
void label_extremal_qqx(
ConstIterator born_begin, ConstIterator born_end,
Iterator first_out
){
// find born quarks
const auto firstquark = std::find_if(
born_begin, born_end-1,
[](Particle const & s){ return (is_anyquark(s)); }
);
assert(firstquark != born_end-1);
const auto secondquark = std::find_if(
firstquark+1, born_end,
[](Particle const & s){ return (is_anyquark(s)); }
);
assert(secondquark != born_end);
assert( ( is_quark(*firstquark) && is_antiquark(*secondquark) )
|| ( is_antiquark(*firstquark) && is_quark(*secondquark) ));
assert(first_out->type == ParticleID::gluon);
assert((first_out+1)->type == ParticleID::gluon);
// copy type from born
first_out->type = firstquark->type;
(first_out+1)->type = secondquark->type;
}
}
void PhaseSpacePoint::label_qqx(Event const & event){
assert(std::is_sorted(begin(outgoing_), end(outgoing_), rapidity_less{}));
assert(filter_partons(outgoing_).size() == outgoing_.size());
if(qqxb_){
label_extremal_qqx( event.outgoing().cbegin(), event.outgoing().cend(),
outgoing_.begin()
);
return;
}
if(qqxf_){ // same as qqxb with reversed order
label_extremal_qqx( event.outgoing().crbegin(), event.outgoing().crend(),
outgoing_.rbegin()
);
return;
}
// central qqx
const auto firstquark = get_first_anyquark_emission(event);
// find jets at FO corresponding to the quarks
// technically this isn't necessary for LO
const auto firstjet_idx = get_back_quark_jet(event, firstquark);
// find corresponding jets after resummation
fastjet::ClusterSequence cs{to_PseudoJet(outgoing_), param_.jet_param.def};
auto const jets = fastjet::sorted_by_rapidity(
cs.inclusive_jets( param_.jet_param.min_pt ));
std::vector<int> const resum_indices{ cs.particle_jet_indices({jets}) };
// assert that jets didn't move
assert(nearby_ep( ( event.jets().cbegin()+firstjet_idx )->rapidity(),
jets[ firstjet_idx ].rapidity(), 1e-2) );
assert(nearby_ep( ( event.jets().cbegin()+firstjet_idx+1 )->rapidity(),
jets[ firstjet_idx+1 ].rapidity(), 1e-2) );
// find last partons in first (central) jet
size_t idx_out = 0;
for(size_t i=resum_indices.size()-2; i>0; --i)
if(resum_indices[i] == firstjet_idx){
idx_out = i;
break;
}
assert(idx_out != 0);
// check that there is sufficient pt in jets from the quarks
const double minpartonjetpt = 1. - param_.max_ext_soft_pt_fraction;
if (outgoing_[idx_out].p.pt()<minpartonjetpt*( event.jets().cbegin()+firstjet_idx )->pt()){
weight_=0.;
status_ = StatusCode::wrong_jets;
return;
}
if (outgoing_[idx_out+1].p.pt()<minpartonjetpt*( event.jets().cbegin()+firstjet_idx+1 )->pt()){
weight_=0.;
status_ = StatusCode::wrong_jets;
return;
}
// check that no additional emission between jets
// such configurations are possible if we have an gluon gets generated
// inside the rapidities of the qqx chain, but clusted to a
// differnet/outside jet. Changing this is non trivial
if(resum_indices[idx_out+1] != resum_indices[idx_out]+1){
weight_=0.;
status_ = StatusCode::gluon_in_qqx;
return;
}
outgoing_[idx_out].type = firstquark->type;
outgoing_[idx_out+1].type = std::next(firstquark)->type;
}
void PhaseSpacePoint::label_quarks(Event const & ev){
const auto WEmit = std::find_if(
begin(ev.outgoing()), end(ev.outgoing()),
[](Particle const & s){ return abs(s.type) == pid::Wp; }
);
if (WEmit != end(ev.outgoing())){
if(!qqxb_) {
const size_t backward_FKL_idx = unob_?1:0;
const auto backward_FKL = std::next(ev.begin_partons(), backward_FKL_idx);
outgoing_[backward_FKL_idx].type = backward_FKL->type;
}
if(!qqxf_) {
const size_t forward_FKL_idx = unof_?1:0;
const auto forward_FKL = std::prev(ev.end_partons(), 1+forward_FKL_idx);
outgoing_.rbegin()[unof_].type = forward_FKL->type;
}
} else {
most_backward_FKL(outgoing_).type = ev.incoming().front().type;
most_forward_FKL(outgoing_).type = ev.incoming().back().type;
}
if(qqxmid_||qqxb_||qqxf_){
label_qqx(ev);
}
}
PhaseSpacePoint::PhaseSpacePoint(
- Event const & ev, PhaseSpacePointConfig conf, HEJ::RNG & ran
+ Event const & ev, PhaseSpacePointConfig conf, RNG & ran
):
unob_{ev.type() == event_type::unob},
unof_{ev.type() == event_type::unof},
qqxb_{ev.type() == event_type::qqxexb},
qqxf_{ev.type() == event_type::qqxexf},
qqxmid_{ev.type() == event_type::qqxmid},
param_{std::move(conf)},
- ran_{ran},
status_{unspecified}
{
weight_ = 1;
const auto & Born_jets = ev.jets();
- const int ng = sample_ng(Born_jets);
+ const int ng = sample_ng(Born_jets, ran);
weight_ /= std::tgamma(ng + 1);
- const int ng_jets = sample_ng_jets(ng, Born_jets);
+ const int ng_jets = sample_ng_jets(ng, Born_jets, ran);
std::vector<fastjet::PseudoJet> out_partons = gen_non_jet(
- ng - ng_jets, CMINPT, param_.jet_param.min_pt
+ ng - ng_jets, CMINPT, param_.jet_param.min_pt, ran
);
int qqxbackjet(-1);
if(qqxmid_){
qqxbackjet = getBackQuarkJet(ev);
}
const auto qperp = std::accumulate(
begin(out_partons), end(out_partons),
fastjet::PseudoJet{}
);
const auto jets = reshuffle(Born_jets, qperp);
if(weight_ == 0.) {
status_ = failed_reshuffle;
return;
}
if(! pass_resummation_cuts(jets)){
status_ = failed_resummation_cuts;
weight_ = 0.;
return;
}
- std::vector<fastjet::PseudoJet> jet_partons = split(jets, ng_jets, qqxbackjet);
+ std::vector<fastjet::PseudoJet> jet_partons = split(
+ jets, ng_jets, qqxbackjet, ran
+ );
if(weight_ == 0.) {
status_ = StatusCode::failed_split;
return;
}
if(qqxmid_){
rescale_qqx_rapidities(
out_partons, jets,
most_backward_FKL(jet_partons).rapidity(),
most_forward_FKL(jet_partons).rapidity(),
qqxbackjet
);
}
else{
rescale_rapidities(
out_partons,
most_backward_FKL(jet_partons).rapidity(),
most_forward_FKL(jet_partons).rapidity()
- );
+ );
}
if(! cluster_jets(out_partons).empty()){
weight_ = 0.;
status_ = StatusCode::empty_jets;
return;
}
std::sort(begin(out_partons), end(out_partons), rapidity_less{});
assert(
std::is_sorted(begin(jet_partons), end(jet_partons), rapidity_less{})
);
const auto first_jet_parton = out_partons.insert(
end(out_partons), begin(jet_partons), end(jet_partons)
);
std::inplace_merge(
begin(out_partons), first_jet_parton, end(out_partons), rapidity_less{}
);
if(! jets_ok(Born_jets, out_partons)){
weight_ = 0.;
status_ = StatusCode::wrong_jets;
return;
}
weight_ *= phase_space_normalisation(Born_jets.size(), out_partons.size());
outgoing_.reserve(out_partons.size() + 1); // one slot for possible A, W, Z, H
for(auto & p: out_partons){
outgoing_.emplace_back(Particle{pid::gluon, std::move(p), {}});
}
assert(!outgoing_.empty());
label_quarks(ev);
if(weight_ == 0.) {
//! @TODO optimise s.t. this is not possible
// status is handled internally
return;
}
copy_AWZH_boson_from(ev);
reconstruct_incoming(ev.incoming());
status_ = StatusCode::good;
}
std::vector<fastjet::PseudoJet> PhaseSpacePoint::gen_non_jet(
- int count, double ptmin, double ptmax
+ int count, double ptmin, double ptmax, RNG & ran
){
// heuristic parameters for pt sampling
const double ptpar = 1.3 + count/5.;
const double temp1 = atan((ptmax - ptmin)/ptpar);
std::vector<fastjet::PseudoJet> partons(count);
for(size_t i = 0; i < (size_t) count; ++i){
- const double r1 = ran_.get().flat();
+ const double r1 = ran.flat();
const double pt = ptmin + ptpar*tan(r1*temp1);
const double temp2 = cos(r1*temp1);
- const double phi = 2*M_PI*ran_.get().flat();
+ const double phi = 2*M_PI*ran.flat();
weight_ *= 2.0*M_PI*pt*ptpar*temp1/(temp2*temp2);
// we don't know the allowed rapidity span yet,
// set a random value to be rescaled later on
- const double y = ran_.get().flat();
+ const double y = ran.flat();
partons[i].reset_PtYPhiM(pt, y, phi);
// Set user index higher than any jet-parton index
// in order to assert that these are not inside jets
partons[i].set_user_index(i + 1 + ng_max);
assert(ptmin-1e-5 <= partons[i].pt() && partons[i].pt() <= ptmax+1e-5);
}
assert(std::all_of(partons.cbegin(), partons.cend(), is_nonjet_parton));
return sorted_by_rapidity(partons);
}
void PhaseSpacePoint::rescale_qqx_rapidities(
std::vector<fastjet::PseudoJet> & out_partons,
std::vector<fastjet::PseudoJet> const & jets,
const double ymin1, const double ymax2,
const int qqxbackjet
){
const double ymax1 = jets[qqxbackjet].rapidity();
const double ymin2 = jets[qqxbackjet+1].rapidity();
constexpr double ep = 1e-7;
const double tot_y = ymax1 - ymin1 + ymax2 - ymin2;
std::vector<std::reference_wrapper<fastjet::PseudoJet>> refpart(
out_partons.begin(), out_partons.end());
double ratio = (ymax1 - ymin1)/tot_y;
const auto gap{ std::find_if(refpart.begin(), refpart.end(),
[ratio](fastjet::PseudoJet p){
return (p.rapidity()>=ratio);} ) };
double ymin = ymin1;
double ymax = ymax1;
double dy = ymax - ymin - 2*ep;
double offset = 0.;
for(auto it_part=refpart.begin(); it_part<refpart.end(); ++it_part){
if(it_part == gap){
ymin = ymin2;
ymax = ymax2;
dy = ymax - ymin - 2*ep;
offset = ratio;
ratio = 1-ratio;
}
fastjet::PseudoJet & part = *it_part;
assert(offset <= part.rapidity() && part.rapidity() < ratio+offset);
const double y = ymin + ep + dy*((part.rapidity()-offset)/ratio);
part.reset_momentum_PtYPhiM(part.pt(), y, part.phi());
weight_ *= tot_y-4.*ep;
assert(ymin <= part.rapidity() && part.rapidity() <= ymax);
}
assert(is_sorted(begin(out_partons), end(out_partons), rapidity_less{}));
}
void PhaseSpacePoint::rescale_rapidities(
std::vector<fastjet::PseudoJet> & partons,
double ymin, double ymax
){
constexpr double ep = 1e-7;
for(auto & parton: partons){
assert(0 <= parton.rapidity() && parton.rapidity() <= 1);
const double dy = ymax - ymin - 2*ep;
const double y = ymin + ep + dy*parton.rapidity();
parton.reset_momentum_PtYPhiM(parton.pt(), y, parton.phi());
weight_ *= dy;
assert(ymin <= parton.rapidity() && parton.rapidity() <= ymax);
}
}
namespace {
template<typename T, typename... Rest>
auto min(T const & a, T const & b, Rest&&... r) {
using std::min;
return min(a, min(b, std::forward<Rest>(r)...));
}
}
double PhaseSpacePoint::probability_in_jet(
std::vector<fastjet::PseudoJet> const & Born_jets
) const{
assert(std::is_sorted(begin(Born_jets), end(Born_jets), rapidity_less{}));
assert(Born_jets.size() >= 2);
const double dy =
Born_jets.back().rapidity() - Born_jets.front().rapidity();
const double R = param_.jet_param.def.R();
const int njets = Born_jets.size();
const double p_J_y_large = (njets-1)*R*R/(2.*dy);
const double p_J_y0 = njets*R/M_PI;
return min(p_J_y_large, p_J_y0, 1.);
}
int PhaseSpacePoint::sample_ng_jets(
- int ng, std::vector<fastjet::PseudoJet> const & Born_jets
+ int ng, std::vector<fastjet::PseudoJet> const & Born_jets, RNG & ran
){
const double p_J = probability_in_jet(Born_jets);
std::binomial_distribution<> bin_dist(ng, p_J);
- const int ng_J = bin_dist(ran_.get());
+ const int ng_J = bin_dist(ran);
weight_ *= std::pow(p_J, -ng_J)*std::pow(1 - p_J, ng_J - ng);
return ng_J;
}
std::vector<fastjet::PseudoJet> PhaseSpacePoint::reshuffle(
std::vector<fastjet::PseudoJet> const & Born_jets,
fastjet::PseudoJet const & q
){
if(q == fastjet::PseudoJet{0, 0, 0, 0}) return Born_jets;
const auto jets = resummation_jet_momenta(Born_jets, q);
if(jets.empty()){
weight_ = 0;
return {};
}
// additional Jacobian to ensure Born integration over delta gives 1
weight_ *= resummation_jet_weight(Born_jets, q);
return jets;
}
std::vector<int> PhaseSpacePoint::distribute_jet_partons(
- int ng_jets, std::vector<fastjet::PseudoJet> const & jets
+ int ng_jets, std::vector<fastjet::PseudoJet> const & jets, RNG & ran
){
size_t first_valid_jet = 0;
size_t num_valid_jets = jets.size();
const double R_eff = 5./3.*param_.jet_param.def.R();
// if there is an unordered jet too far away from the FKL jets
// then extra gluon constituents of the unordered jet would
// violate the FKL rapidity ordering
if((unob_||qqxb_) && jets[0].delta_R(jets[1]) > R_eff){
++first_valid_jet;
--num_valid_jets;
}
else if((unof_||qqxf_) && jets[jets.size()-1].delta_R(jets[jets.size()-2]) > R_eff){
--num_valid_jets;
}
std::vector<int> np(jets.size(), 1);
for(int i = 0; i < ng_jets; ++i){
- ++np[first_valid_jet + ran_.get().flat() * num_valid_jets];
+ ++np[first_valid_jet + ran.flat() * num_valid_jets];
}
weight_ *= std::pow(num_valid_jets, ng_jets);
return np;
}
#ifndef NDEBUG
namespace{
bool tagged_FKL_backward(
std::vector<fastjet::PseudoJet> const & jet_partons
){
return std::find_if(
begin(jet_partons), end(jet_partons),
[](fastjet::PseudoJet const & p){
return p.user_index() == backward_FKL_idx;
}
) != end(jet_partons);
}
bool tagged_FKL_forward(
std::vector<fastjet::PseudoJet> const & jet_partons
){
// the most forward FKL parton is most likely near the end of jet_partons;
// start search from there
return std::find_if(
jet_partons.rbegin(), jet_partons.rend(),
[](fastjet::PseudoJet const & p){
return p.user_index() == forward_FKL_idx;
}
) != jet_partons.rend();
}
bool tagged_FKL_extremal(
std::vector<fastjet::PseudoJet> const & jet_partons
){
return tagged_FKL_backward(jet_partons) && tagged_FKL_forward(jet_partons);
}
} // namespace anonymous
#endif
std::vector<fastjet::PseudoJet> PhaseSpacePoint::split(
std::vector<fastjet::PseudoJet> const & jets,
- int ng_jets,
- size_t qqxbackjet
+ int ng_jets, size_t qqxbackjet, RNG & ran
){
- return split(jets, distribute_jet_partons(ng_jets, jets), qqxbackjet);
+ return split(
+ jets, distribute_jet_partons(ng_jets, jets, ran), qqxbackjet, ran);
}
bool PhaseSpacePoint::pass_extremal_cuts(
fastjet::PseudoJet const & ext_parton,
fastjet::PseudoJet const & jet
) const{
if(ext_parton.pt() < param_.min_extparton_pt) return false;
return (ext_parton - jet).pt()/jet.pt() < param_.max_ext_soft_pt_fraction;
}
std::vector<fastjet::PseudoJet> PhaseSpacePoint::split(
std::vector<fastjet::PseudoJet> const & jets,
std::vector<int> const & np,
- size_t qqxbackjet
+ size_t qqxbackjet,
+ RNG & ran
){
assert(! jets.empty());
assert(jets.size() == np.size());
assert(pass_resummation_cuts(jets));
const size_t most_backward_FKL_idx = 0 + unob_ + qqxb_;
const size_t most_forward_FKL_idx = jets.size() - 1 - unof_ - qqxf_;
const auto & jet = param_.jet_param;
- const JetSplitter jet_splitter{jet.def, jet.min_pt, ran_};
+ const JetSplitter jet_splitter{jet.def, jet.min_pt};
std::vector<fastjet::PseudoJet> jet_partons;
// randomly distribute jet gluons among jets
for(size_t i = 0; i < jets.size(); ++i){
- auto split_res = jet_splitter.split(jets[i], np[i]);
+ auto split_res = jet_splitter.split(jets[i], np[i], ran);
weight_ *= split_res.weight;
if(weight_ == 0) return {};
assert(
std::all_of(
begin(split_res.constituents), end(split_res.constituents),
is_jet_parton
)
);
const auto first_new_parton = jet_partons.insert(
end(jet_partons),
begin(split_res.constituents), end(split_res.constituents)
);
// mark uno and extremal FKL emissions here so we can check
// their position once all emissions are generated
// also mark qqxmid partons, and apply appropriate pt cut.
auto extremal = end(jet_partons);
if (i == most_backward_FKL_idx){ //FKL backward emission
extremal = std::min_element(
first_new_parton, end(jet_partons), rapidity_less{}
);
extremal->set_user_index(backward_FKL_idx);
}
else if(((unob_ || qqxb_) && i == 0)){
// unordered/qqxb
extremal = std::min_element(
first_new_parton, end(jet_partons), rapidity_less{}
);
extremal->set_user_index((unob_)?unob_idx:qqxb_idx);
}
else if (i == most_forward_FKL_idx){
extremal = std::max_element(
first_new_parton, end(jet_partons), rapidity_less{}
);
extremal->set_user_index(forward_FKL_idx);
}
else if(((unof_ || qqxf_) && i == jets.size() - 1)){
// unordered/qqxf
extremal = std::max_element(
first_new_parton, end(jet_partons), rapidity_less{}
);
extremal->set_user_index((unof_)?unof_idx:qqxf_idx);
}
else if((qqxmid_ && i == qqxbackjet)){
extremal = std::max_element(
first_new_parton, end(jet_partons), rapidity_less{}
);
extremal->set_user_index(qqxmid1_idx);
}
else if((qqxmid_ && i == qqxbackjet+1)){
extremal = std::min_element(
first_new_parton, end(jet_partons), rapidity_less{}
);
extremal->set_user_index(qqxmid2_idx);
}
if(
extremal != end(jet_partons)
&& !pass_extremal_cuts(*extremal, jets[i])
){
weight_ = 0;
return {};
}
}
assert(tagged_FKL_extremal(jet_partons));
std::sort(begin(jet_partons), end(jet_partons), rapidity_less{});
if(
!extremal_ok(jet_partons)
|| !split_preserved_jets(jets, jet_partons)
){
weight_ = 0.;
return {};
}
return jet_partons;
}
bool PhaseSpacePoint::extremal_ok(
std::vector<fastjet::PseudoJet> const & partons
) const{
assert(std::is_sorted(begin(partons), end(partons), rapidity_less{}));
if(unob_ && partons.front().user_index() != unob_idx) return false;
if(unof_ && partons.back().user_index() != unof_idx) return false;
if(qqxb_ && partons.front().user_index() != qqxb_idx) return false;
if(qqxf_ && partons.back().user_index() != qqxf_idx) return false;
return
most_backward_FKL(partons).user_index() == backward_FKL_idx
&& most_forward_FKL(partons).user_index() == forward_FKL_idx;
}
bool PhaseSpacePoint::split_preserved_jets(
std::vector<fastjet::PseudoJet> const & jets,
std::vector<fastjet::PseudoJet> const & jet_partons
) const{
assert(std::is_sorted(begin(jets), end(jets), rapidity_less{}));
const auto split_jets = cluster_jets(jet_partons);
// this can happen if two overlapping jets
// are both split into more than one parton
if(split_jets.size() != jets.size()) return false;
for(size_t i = 0; i < split_jets.size(); ++i){
// this can happen if there are two overlapping jets
// and a parton is assigned to the "wrong" jet
if(!nearby_ep(jets[i].rapidity(), split_jets[i].rapidity(), 1e-2)){
return false;
}
}
return true;
}
template<class Particle>
Particle const & PhaseSpacePoint::most_backward_FKL(
std::vector<Particle> const & partons
) const{
return partons[0 + unob_ + qqxb_];
}
template<class Particle>
Particle const & PhaseSpacePoint::most_forward_FKL(
std::vector<Particle> const & partons
) const{
const size_t idx = partons.size() - 1 - unof_ - qqxf_;
assert(idx < partons.size());
return partons[idx];
}
template<class Particle>
Particle & PhaseSpacePoint::most_backward_FKL(
std::vector<Particle> & partons
) const{
return partons[0 + unob_ + qqxb_];
}
template<class Particle>
Particle & PhaseSpacePoint::most_forward_FKL(
std::vector<Particle> & partons
) const{
const size_t idx = partons.size() - 1 - unof_ - qqxf_;
assert(idx < partons.size());
return partons[idx];
}
bool PhaseSpacePoint::contains_idx(
fastjet::PseudoJet const & jet, fastjet::PseudoJet const & parton
) const {
auto const & constituents = jet.constituents();
const int idx = parton.user_index();
const bool injet = std::find_if(
begin(constituents), end(constituents),
[idx](fastjet::PseudoJet const & con){return con.user_index() == idx;}
) != end(constituents);
const double minpartonjetpt = 1. - param_.max_ext_soft_pt_fraction;
return ((parton.pt()>minpartonjetpt*jet.pt())&&injet);
}
bool PhaseSpacePoint::jets_ok(
std::vector<fastjet::PseudoJet> const & Born_jets,
std::vector<fastjet::PseudoJet> const & partons
) const{
fastjet::ClusterSequence cs(partons, param_.jet_param.def);
const auto jets = sorted_by_rapidity(cs.inclusive_jets(param_.jet_param.min_pt));
if(jets.size() != Born_jets.size()) return false;
int in_jet = 0;
for(size_t i = 0; i < jets.size(); ++i){
assert(jets[i].has_constituents());
for(auto && parton: jets[i].constituents()){
if(is_nonjet_parton(parton)) return false;
}
in_jet += jets[i].constituents().size();
}
const int expect_in_jet = std::count_if(
partons.cbegin(), partons.cend(), is_jet_parton
);
if(in_jet != expect_in_jet) return false;
// note that PseudoJet::contains does not work here
if(! (
contains_idx(most_backward_FKL(jets), most_backward_FKL(partons))
&& contains_idx(most_forward_FKL(jets), most_forward_FKL(partons))
)) return false;
if(unob_ && !contains_idx(jets.front(), partons.front())) return false;
if(qqxb_ && !contains_idx(jets.front(), partons.front())) return false;
if(unof_ && !contains_idx(jets.back(), partons.back())) return false;
if(qqxf_ && !contains_idx(jets.back(), partons.back())) return false;
for(size_t i = 0; i < jets.size(); ++i){
assert(nearby_ep(jets[i].rapidity(), Born_jets[i].rapidity(), 1e-2));
}
return true;
}
void PhaseSpacePoint::reconstruct_incoming(
std::array<Particle, 2> const & Born_incoming
){
std::tie(incoming_[0].p, incoming_[1].p) = incoming_momenta(outgoing_);
for(size_t i = 0; i < incoming_.size(); ++i){
incoming_[i].type = Born_incoming[i].type;
}
assert(momentum_conserved());
}
double PhaseSpacePoint::phase_space_normalisation(
int num_Born_jets, int num_out_partons
) const{
return pow(16*pow(M_PI,3), num_Born_jets - num_out_partons);
}
bool PhaseSpacePoint::momentum_conserved() const{
fastjet::PseudoJet diff;
for(auto const & in: incoming()) diff += in.p;
const double norm = diff.E();
for(auto const & out: outgoing()) diff -= out.p;
return nearby(diff, fastjet::PseudoJet{}, norm);
}
} //namespace HEJ
diff --git a/src/bin/HEJ.cc b/src/bin/HEJ.cc
index ccc20ff..37067c1 100644
--- a/src/bin/HEJ.cc
+++ b/src/bin/HEJ.cc
@@ -1,375 +1,376 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include <array>
#include <chrono>
#include <iostream>
#include <limits>
#include <memory>
#include <numeric>
#include "yaml-cpp/yaml.h"
#include "fastjet/ClusterSequence.hh"
#include "HEJ/CombinedEventWriter.hh"
#include "HEJ/Config.hh"
#include "HEJ/CrossSectionAccumulator.hh"
#include "HEJ/Event.hh"
#include "HEJ/EventReader.hh"
#include "HEJ/BufferedEventReader.hh"
#include "HEJ/EventReweighter.hh"
#include "HEJ/get_analysis.hh"
#include "HEJ/make_RNG.hh"
#include "HEJ/optional.hh"
#include "HEJ/ProgressBar.hh"
#include "HEJ/stream.hh"
#include "HEJ/Unweighter.hh"
#include "HEJ/Version.hh"
#include "HEJ/YAMLreader.hh"
HEJ::Config load_config(char const * filename){
try{
return HEJ::load_config(filename);
}
catch(std::exception const & exc){
std::cerr << "Error: " << exc.what() << '\n';
std::exit(EXIT_FAILURE);
}
}
std::unique_ptr<HEJ::Analysis> get_analysis(
YAML::Node const & parameters, LHEF::HEPRUP const & heprup
){
try{
return HEJ::get_analysis(parameters, heprup);
}
catch(std::exception const & exc){
std::cerr << "Failed to load analysis: " << exc.what() << '\n';
std::exit(EXIT_FAILURE);
}
}
// unique_ptr is a workaround:
// HEJ::optional is a better fit, but gives spurious errors with g++ 7.3.0
std::unique_ptr<HEJ::ProgressBar<double>> make_progress_bar(
std::vector<double> const & xs
) {
if(xs.empty()) return {};
const double Born_xs = std::accumulate(begin(xs), end(xs), 0.);
return std::make_unique<HEJ::ProgressBar<double>>(std::cout, Born_xs);
}
std::string time_to_string(const time_t time){
char s[30];
struct tm * p = localtime(&time);
strftime(s, 30, "%a %b %d %Y %H:%M:%S", p);
return s;
}
HEJ::Event to_event(
LHEF::HEPEUP const & hepeup,
HEJ::JetParameters const & fixed_order_jets
) {
HEJ::Event::EventData event_data{hepeup};
event_data.reconstruct_intermediate();
return HEJ::Event{
std::move(event_data).cluster(
fixed_order_jets.def, fixed_order_jets.min_pt
)
};
}
void unweight(
HEJ::Unweighter & unweighter,
HEJ::WeightType weight_type,
std::vector<HEJ::Event> & events,
HEJ::RNG & ran
) {
if(weight_type == HEJ::WeightType::unweighted_resum){
unweighter.set_cut_to_maxwt(events);
}
events.erase(
unweighter.unweight(begin(events), end(events), ran),
end(events)
);
}
// peek up to nevents events from reader
std::vector<LHEF::HEPEUP> peek_events(
HEJ::BufferedEventReader & reader,
const int nevents
) {
std::vector<LHEF::HEPEUP> events;
while(
static_cast<int>(events.size()) < nevents
&& reader.read_event()
) {
events.emplace_back(reader.hepeup());
}
// put everything back into the reader
for(auto it = rbegin(events); it != rend(events); ++it) {
reader.emplace(*it);
}
return events;
}
void append_resummed_events(
std::vector<HEJ::Event> & resummation_events,
HEJ::EventReweighter & reweighter,
LHEF::HEPEUP const & hepeup,
const int trials,
HEJ::JetParameters const & fixed_order_jets
) {
const HEJ::Event FO_event = to_event(hepeup, fixed_order_jets);
if(reweighter.treatment(FO_event.type()) != HEJ::EventTreatment::reweight) {
return;
}
const auto resummed = reweighter.reweight(FO_event, trials);
resummation_events.insert(
end(resummation_events),
begin(resummed), end(resummed)
);
}
void train(
HEJ::Unweighter & unweighter,
HEJ::BufferedEventReader & reader,
HEJ::EventReweighter & reweighter,
const int total_trials,
const double max_dev,
double reweight_factor,
HEJ::JetParameters const & fixed_order_jets
) {
std::cout << "Reading up to " << total_trials << " training events...\n";
auto FO_events = peek_events(reader, total_trials);
if(FO_events.empty()) {
throw std::runtime_error{
"No events generated to calibrate the unweighting weight!"
"Please increase the number \"trials\" or deactivate the unweighting."
};
}
const int trials = total_trials/FO_events.size();
// adjust reweight factor so that the overall normalisation
// is the same as in the full run
reweight_factor *= trials;
for(auto & hepeup: FO_events) {
hepeup.XWGTUP *= reweight_factor;
}
std::cout << "Training unweighter with "
<< trials << '*' << FO_events.size() << " events\n";
auto progress = HEJ::ProgressBar<int>{
std::cout, static_cast<int>(FO_events.size())
};
std::vector<HEJ::Event> resummation_events;
for(auto const & hepeup: FO_events) {
append_resummed_events(
resummation_events,
reweighter, hepeup, trials, fixed_order_jets
);
++progress;
}
unweighter.set_cut_to_peakwt(resummation_events, max_dev);
std::cout << "\nUnweighting events with weight up to "
<< unweighter.get_cut() << '\n';
}
int main(int argn, char** argv) {
using clock = std::chrono::system_clock;
if (argn != 3) {
std::cerr << "\n# Usage:\n."<< argv[0] <<" config_file input_file\n\n";
return EXIT_FAILURE;
}
const auto start_time = clock::now();
{
std::cout << "Starting " << HEJ::Version::package_name_full()
<< ", revision " << HEJ::Version::revision() << " ("
<< time_to_string(clock::to_time_t(start_time)) << ")" << std::endl;
}
fastjet::ClusterSequence::print_banner();
// read configuration
const HEJ::Config config = load_config(argv[1]);
auto reader = HEJ::make_reader(argv[2]);
assert(reader);
auto heprup{ reader->heprup() };
heprup.generators.emplace_back(LHEF::XMLTag{});
heprup.generators.back().name = HEJ::Version::package_name();
heprup.generators.back().version = HEJ::Version::String();
std::unique_ptr<HEJ::Analysis> analysis = get_analysis(
config.analysis_parameters, heprup
);
assert(analysis != nullptr);
HEJ::CombinedEventWriter writer{config.output, std::move(heprup)};
double global_reweight = 1.;
const auto & max_events = config.max_events;
// if we need the event number:
if(std::abs(heprup.IDWTUP) == 4 || std::abs(heprup.IDWTUP) == 1 || max_events){
// try to read from LHE head
auto input_events{reader->number_events()};
if(!input_events) {
// else count manually
auto t_reader = HEJ::make_reader(argv[2]);
input_events = 0;
while(t_reader->read_event()) ++(*input_events);
}
if(std::abs(heprup.IDWTUP) == 4 || std::abs(heprup.IDWTUP) == 1){
// IDWTUP 4 or 1 assume average(weight)=xs, but we need sum(weights)=xs
std::cout << "Found IDWTUP " << heprup.IDWTUP << ": "
<< "assuming \"cross section = average weight\".\n"
<< "converting to \"cross section = sum of weights\" ";
global_reweight /= *input_events;
}
if(max_events && (*input_events > *max_events)){
// maximal number of events given
global_reweight *= *input_events/static_cast<double>(*max_events);
std::cout << "Processing " << *max_events
<< " out of " << *input_events << " events\n";
}
}
HEJ::ScaleGenerator scale_gen{
config.scales.base,
config.scales.factors,
config.scales.max_ratio
};
- auto ran = HEJ::make_RNG(config.rng.name, config.rng.seed);
+ std::shared_ptr<HEJ::RNG> ran{
+ HEJ::make_RNG(config.rng.name, config.rng.seed)};
assert(ran != nullptr);
HEJ::EventReweighter hej{
reader->heprup(),
std::move(scale_gen),
to_EventReweighterConfig(config),
- *ran
+ ran
};
HEJ::optional<HEJ::Unweighter> unweighter{};
if(config.weight_type != HEJ::WeightType::weighted) {
unweighter = HEJ::Unweighter{};
}
if(config.weight_type == HEJ::WeightType::partially_unweighted) {
HEJ::BufferedEventReader buffered_reader{std::move(reader)};
assert(config.unweight_config);
train(
*unweighter,
buffered_reader,
hej,
config.unweight_config->trials,
config.unweight_config->max_dev,
global_reweight/config.trials,
config.fixed_order_jets
);
reader = std::make_unique<HEJ::BufferedEventReader>(
std::move(buffered_reader)
);
}
// status infos & eye candy
size_t nevent = 0;
std::array<int, HEJ::event_type::last_type + 1>
nevent_type{0}, nfailed_type{0};
auto progress = make_progress_bar(reader->heprup().XSECUP);
HEJ::CrossSectionAccumulator xs;
std::map<HEJ::StatusCode, int> status_counter;
size_t total_trials = 0;
size_t total_resum = 0;
// Loop over the events in the input file
while(reader->read_event() && (!max_events || nevent < *max_events) ){
++nevent;
// reweight events so that the total cross section is conserved
auto hepeup = reader->hepeup();
hepeup.XWGTUP *= global_reweight;
const auto FO_event = to_event(hepeup, config.fixed_order_jets);
if(FO_event.central().weight == 0) {
static const bool warned_once = [argv,nevent](){
std::cerr
<< "WARNING: event number " << nevent
<< " in " << argv[2] << " has zero weight. "
"Ignoring this and all further events with vanishing weight.\n";
return true;
}();
(void) warned_once; // shut up compiler warnings
continue;
}
auto resummed_events{ hej.reweight(FO_event, config.trials) };
// some bookkeeping
for(auto const & s: hej.status())
++status_counter[s];
total_trials+=hej.status().size();
++nevent_type[FO_event.type()];
if(resummed_events.empty()) ++nfailed_type[FO_event.type()];
if(unweighter) {
unweight(*unweighter, config.weight_type, resummed_events, *ran);
}
// analysis
for(auto & ev: resummed_events){
//TODO: move pass_cuts to after phase space point generation
if(analysis->pass_cuts(ev, FO_event)){
analysis->fill(ev, FO_event);
writer.write(ev);
} else {
ev.parameters()*=0; // do not use discarded events afterwards
}
}
xs.fill_correlated(resummed_events);
total_resum += resummed_events.size();
if(progress) progress->increment(FO_event.central().weight);
} // main event loop
std::cout << '\n';
analysis->finalise();
using namespace HEJ::event_type;
std::cout<< "Events processed: " << nevent << " (" << total_resum << " resummed)"<< '\n';
std::cout << '\t' << name(EventType::first_type) << ": "
<< nevent_type[EventType::first_type]
<< ", failed to reconstruct " << nfailed_type[EventType::first_type]
<< '\n';
for(auto i=EventType::first_type+1; i<=EventType::last_type; i*=2){
std::cout << '\t' << name(static_cast<EventType>(i)) << ": "
<< nevent_type[i]
<< ", failed to reconstruct " << nfailed_type[i]
<< '\n';
}
std::cout << '\n' << xs << '\n';
std::cout << "Generation statistic: "
<< status_counter[HEJ::StatusCode::good] << "/" << total_trials
<< " trials successful.\n";
for(auto && entry: status_counter){
const double fraction = static_cast<double>(entry.second)/total_trials;
const int percent = std::round(100*fraction);
std::cout << std::left << std::setw(17) << (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 << "] " <<std::setw(2)<<std::right<< percent << "%\n";
}
std::chrono::duration<double> run_time = (clock::now() - start_time);
std::cout << "\nFinished " << HEJ::Version::package_name() << " at "
<< time_to_string(clock::to_time_t(clock::now()))
<< "\n=> Runtime: " << run_time.count() << " sec ("
<< nevent/run_time.count() << " Events/sec).\n";
return EXIT_SUCCESS;
}
diff --git a/t/check_res.cc b/t/check_res.cc
index 5b6cd59..3e74e19 100644
--- a/t/check_res.cc
+++ b/t/check_res.cc
@@ -1,168 +1,168 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include <cmath>
#include <iostream>
#include "LHEF/LHEF.h"
#include "HEJ/CrossSectionAccumulator.hh"
#include "HEJ/Event.hh"
#include "HEJ/EventReweighter.hh"
#include "HEJ/Mixmax.hh"
#include "HEJ/stream.hh"
#include "hej_test.hh"
namespace{
const fastjet::JetDefinition jet_def{fastjet::kt_algorithm, 0.4};
const fastjet::JetDefinition Born_jet_def{jet_def};
constexpr double Born_jetptmin = 30;
constexpr double max_ext_soft_pt_fraction = 0.1;
constexpr double jetptmin = 35;
constexpr bool log_corr = false;
const HEJ::ParticleProperties Wprop{80.385, 2.085};
const HEJ::ParticleProperties Zprop{91.187, 2.495};
const HEJ::ParticleProperties Hprop{125, 0.004165};
constexpr double vev = 246.2196508;
using EventTreatment = HEJ::EventTreatment;
using namespace HEJ::event_type;
HEJ::EventTreatMap treat{
{no_2_jets, EventTreatment::discard},
{bad_final_state, EventTreatment::discard},
{non_resummable, EventTreatment::discard},
{unof, EventTreatment::discard},
{unob, EventTreatment::discard},
{qqxexb, EventTreatment::discard},
{qqxexf, EventTreatment::discard},
{qqxmid, EventTreatment::discard},
{FKL, EventTreatment::reweight}
};
/// true if colour is allowed for particle
bool correct_colour(HEJ::Particle const & part){
if(HEJ::is_AWZH_boson(part) && !part.colour) return true;
if(!part.colour) return false;
int const colour = part.colour->first;
int const anti_colour = part.colour->second;
if(part.type == HEJ::ParticleID::gluon)
return colour != anti_colour && colour > 0 && anti_colour > 0;
if(HEJ::is_quark(part))
return anti_colour == 0 && colour > 0;
return colour == 0 && anti_colour > 0;
}
bool correct_colour(HEJ::Event const & ev){
if(!HEJ::event_type::is_resummable(ev.type()))
return true;
for(auto const & part: ev.incoming()){
if(!correct_colour(part))
return false;
}
for(auto const & part: ev.outgoing()){
if(!correct_colour(part))
return false;
}
return true;
}
}
int main(int argn, char** argv) {
if(argn == 5 && std::string(argv[4]) == "unof"){
--argn;
treat[unof] = EventTreatment::reweight;
treat[unob] = EventTreatment::discard;
treat[FKL] = EventTreatment::discard;
}
if(argn == 5 && std::string(argv[4]) == "unob"){
--argn;
treat[unof] = EventTreatment::discard;
treat[unob] = EventTreatment::reweight;
treat[FKL] = EventTreatment::discard;
}
else if(argn == 5 && std::string(argv[4]) == "splitf"){
--argn;
treat[qqxexb] = EventTreatment::discard;
treat[qqxexf] = EventTreatment::reweight;
treat[FKL] = EventTreatment::discard;
}
else if(argn == 5 && std::string(argv[4]) == "splitb"){
--argn;
treat[qqxexb] = EventTreatment::reweight;
treat[qqxexf] = EventTreatment::discard;
treat[FKL] = EventTreatment::discard;
}
else if(argn == 5 && std::string(argv[4]) == "qqxmid"){
--argn;
treat[qqxmid] = EventTreatment::reweight;
treat[FKL] = EventTreatment::discard;
}
if(argn != 4){
std::cerr << "Usage: check_res eventfile xsection tolerance [uno]";
return EXIT_FAILURE;
}
const double xsec_ref = std::stod(argv[2]);
const double tolerance = std::stod(argv[3]);
HEJ::istream in{argv[1]};
LHEF::Reader reader{in};
HEJ::PhaseSpacePointConfig psp_conf;
psp_conf.jet_param = HEJ::JetParameters{jet_def, jetptmin};
psp_conf.max_ext_soft_pt_fraction = max_ext_soft_pt_fraction;
HEJ::MatrixElementConfig ME_conf;
ME_conf.log_correction = log_corr;
ME_conf.Higgs_coupling = HEJ::HiggsCouplingSettings{};
ME_conf.ew_parameters.set_vevWZH(vev, Wprop, Zprop, Hprop);
HEJ::EventReweighterConfig conf;
conf.psp_config = std::move(psp_conf);
conf.ME_config = std::move(ME_conf);
conf.treat = treat;
reader.readEvent();
const bool has_Higgs = std::find(
begin(reader.hepeup.IDUP),
end(reader.hepeup.IDUP),
25
) != end(reader.hepeup.IDUP);
const double mu = has_Higgs?125.:91.188;
HEJ::ScaleGenerator scale_gen{
{{std::to_string(mu), HEJ::FixedScale{mu}}}, {}, 1.
};
- HEJ::Mixmax ran{};
+ std::shared_ptr<HEJ::RNG> ran{std::make_shared<HEJ::Mixmax>()};
HEJ::EventReweighter hej{reader.heprup, std::move(scale_gen), conf, ran};
HEJ::CrossSectionAccumulator xs;
do{
auto ev_data = HEJ::Event::EventData{reader.hepeup};
shuffle_particles(ev_data);
ev_data.reconstruct_intermediate();
HEJ::Event ev{
ev_data.cluster(
Born_jet_def, Born_jetptmin
)
};
auto resummed_events = hej.reweight(ev, 100);
for(auto const & ev: resummed_events) {
ASSERT(correct_colour(ev));
ASSERT(std::isfinite(ev.central().weight));
// we fill the xs uncorrelated since we only want to test the uncertainty
// of the resummation
xs.fill(ev);
}
} while(reader.readEvent());
const double xsec = xs.total().value;
const double xsec_err = std::sqrt(xs.total().error);
const double significance =
std::abs(xsec - xsec_ref) / std::sqrt( xsec_err*xsec_err + tolerance*tolerance );
std::cout << xsec_ref << " +/- " << tolerance << " ~ "
<< xsec << " +- " << xsec_err << " => " << significance << " sigma\n";
if(significance > 3.){
std::cerr << "Cross section is off by over 3 sigma!\n";
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
diff --git a/t/test_scale_arithmetics.cc b/t/test_scale_arithmetics.cc
index 3d6dcfa..268fc9e 100644
--- a/t/test_scale_arithmetics.cc
+++ b/t/test_scale_arithmetics.cc
@@ -1,94 +1,95 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include <fstream>
#include <algorithm>
#include "LHEF/LHEF.h"
#include "HEJ/EventReweighter.hh"
#include "HEJ/make_RNG.hh"
#include "HEJ/Event.hh"
#include "HEJ/YAMLreader.hh"
#include "HEJ/stream.hh"
#include "hej_test.hh"
constexpr double ep = 1e-13;
void dump(HEJ::Event const & ev){
{
LHEF::Writer writer{std::cout};
std::cout << std::setprecision(6);
writer.hepeup = to_HEPEUP(std::move(ev), nullptr);
writer.writeEvent();
}
std::cout << "Rapidity ordering:\n";
for(const auto & part: ev.outgoing()){
std::cout << std::setw(2) << part.type << ": "<< std::setw(7) << part.rapidity() << std::endl;
}
}
int main(int argn, char** argv){
if(argn != 3){
std::cerr << "\n# Usage:\n."<< argv[0] <<" config.yml input_file.lhe\n\n";
return EXIT_FAILURE;
}
HEJ::Config config = HEJ::load_config(argv[1]);
config.scales = HEJ::to_ScaleConfig(
YAML::Load("scales: [H_T, 1 * H_T, 2/2 * H_T, 2*H_T/2, H_T/2*2, H_T/2/2*4, H_T*H_T/H_T]")
);
HEJ::istream in{argv[2]};
LHEF::Reader reader{in};
- auto ran = HEJ::make_RNG(config.rng.name, config.rng.seed);
+ std::shared_ptr<HEJ::RNG> ran{
+ HEJ::make_RNG(config.rng.name, config.rng.seed)};
HEJ::ScaleGenerator scale_gen{
config.scales.base,
config.scales.factors,
config.scales.max_ratio
};
HEJ::EventReweighter resum{
reader.heprup,
std::move(scale_gen),
to_EventReweighterConfig(config),
- *ran
+ ran
};
size_t i = 0;
while(reader.readEvent()){
++i;
HEJ::Event::EventData data{reader.hepeup};
shuffle_particles(data);
HEJ::Event event{
data.cluster(
config.resummation_jets.def,
config.resummation_jets.min_pt
)
};
auto resummed = resum.reweight(event, config.trials);
for(auto && ev: resummed) {
for(auto &&var: ev.variations()) {
if(std::abs(var.muf - ev.central().muf) > ep) {
std::cerr
<< std::setprecision(15)
<< "unequal scales: " << var.muf
<< " != " << ev.central().muf << '\n'
<< "in resummed event:\n";
dump(ev);
std::cerr << "\noriginal event:\n";
dump(event);
return EXIT_FAILURE;
}
}
}
}
return EXIT_SUCCESS;
}

File Metadata

Mime Type
text/x-diff
Expires
Sun, Feb 23, 2:59 PM (52 m, 33 s)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4485158
Default Alt Text
(121 KB)

Event Timeline