Page MenuHomeHEPForge
Files F8724285

No OneTemporary
Actions

View Options

diff --git a/Changes.md b/Changes.md
index d1f98d3..aea4b72 100644
--- a/Changes.md
+++ b/Changes.md
@@ -1,31 +1,33 @@
# Version 2.0
## 2.X.0
* Made `MatrixElement::tree_kin` and `MatrixElement::tree_param` public
* Allow multiplication and division of multiple scale functions e.g.
`H_T/2*m_j1j2`
* Follow HepMC convention for particle Status codes: incoming = 11,
decaying = 2, outgoing = 1 (unchanged)
+* New class `CrossSectionAccumulator` to keep track of Cross Section of the
+ different subproccess
## 2.0.4
* Fixed wrong path of `HEJ_INCLUDE_DIR` in `hej-config.cmake`
* Correctly include rivet headers
* New `EXCLUDE_package` variable in `cmake` to not interface to specific
packages
* Consistent search and include for side packages in `cmake`
## 2.0.3
* Fixed parsing of (numerical factor) * (base scale) in configuration
* Don't change scale names, but sanitise Rivet output file names instead
## 2.0.2
* Changed scale names to `"_over_"` and `"_times_"` for proper file names (was
`"/"` and `"*"` before)
## 2.0.1
* Fixed name of fixed-order generator in error message.
diff --git a/FixedOrderGen/include/CrossSectionAccumulator.hh b/FixedOrderGen/include/CrossSectionAccumulator.hh
deleted file mode 100644
index 00f20c6..0000000
--- a/FixedOrderGen/include/CrossSectionAccumulator.hh
+++ /dev/null
@@ -1,43 +0,0 @@
-#pragma once
-
-#include <map>
-
-#include "HEJ/Event.hh"
-#include "HEJ/event_types.hh"
-
-namespace HEJFOG {
- template<typename T>
- struct WithError {
- T value = T{};
- T error = T{};
- };
-
- class CrossSectionAccumulator {
- public:
- void fill(HEJ::Event const & ev) {
- const double wt = ev.central().weight;
- auto & entry = xs_[ev.type()];
- entry.value += wt;
- entry.error += wt*wt;
- total_.value += wt;
- total_.error += wt*wt;
- }
-
- auto begin() const {
- return std::begin(xs_);
- }
-
- auto end() const {
- return std::end(xs_);
- }
-
- WithError<double> total() const {
- return total_;
- }
-
- private:
- std::map<HEJ::event_type::EventType, WithError<double>> xs_;
- WithError<double> total_;
- };
-
-}
diff --git a/FixedOrderGen/src/main.cc b/FixedOrderGen/src/main.cc
index 7e0cdf3..6d55edb 100644
--- a/FixedOrderGen/src/main.cc
+++ b/FixedOrderGen/src/main.cc
@@ -1,246 +1,233 @@
/**
* Name: main.cc
* Authors: Jeppe R. Andersen
*/
#include <algorithm>
#include <chrono>
#include <fstream>
#include <iostream>
#include <map>
#include <memory>
#include "yaml-cpp/yaml.h"
#include "LHEF/LHEF.h"
#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 "config.hh"
-#include "CrossSectionAccumulator.hh"
#include "EventGenerator.hh"
#include "PhaseSpacePoint.hh"
#include "Unweighter.hh"
#include "Version.hh"
namespace{
constexpr auto banner =
" __ ___ __ ______ __ __ \n"
" / / / (_)___ _/ /_ / ____/___ ___ _________ ___ __ / /__ / /______ \n"
" / /_/ / / __ `/ __ \\ / __/ / __ \\/ _ \\/ ___/ __ `/ / / / __ / / _ \\/ __/ ___/ \n"
" / __ / / /_/ / / / / / /___/ / / / __/ / / /_/ / /_/ / / /_/ / __/ /_(__ ) \n"
" /_/ /_/_/\\__, /_/ /_/ /_____/_/ /_/\\___/_/ \\__, /\\__, / \\____/\\___/\\__/____/ \n"
" ____///__/ __ ____ ///__//____/ ______ __ \n"
" / ____(_) _____ ____/ / / __ \\_________/ /__ _____ / ____/__ ____ ___ _________ _/ /_____ _____\n"
" / /_ / / |/_/ _ \\/ __ / / / / / ___/ __ / _ \\/ ___/ / / __/ _ \\/ __ \\/ _ \\/ ___/ __ `/ __/ __ \\/ ___/\n"
" / __/ / /> </ __/ /_/ / / /_/ / / / /_/ / __/ / / /_/ / __/ / / / __/ / / /_/ / /_/ /_/ / / \n"
" /_/ /_/_/|_|\\___/\\__,_/ \\____/_/ \\__,_/\\___/_/ \\____/\\___/_/ /_/\\___/_/ \\__,_/\\__/\\____/_/ \n"
;
constexpr double invGeV2_to_pb = 389379292.;
constexpr long long max_warmup_events = 10000;
}
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
){
try{
return HEJ::get_analysis(parameters);
}
catch(std::exception const & exc){
std::cerr << "Failed to load analysis: " << exc.what() << '\n';
std::exit(EXIT_FAILURE);
}
}
void rescale_weights(HEJ::Event & ev, double factor) {
ev.central().weight *= factor;
for(auto & var: ev.variations()){
var.weight *= factor;
}
}
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]);
std::unique_ptr<HEJ::Analysis> analysis = get_analysis(
config.analysis_parameters
);
assert(analysis != nullptr);
auto 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.particles_properties,
config.Higgs_coupling,
*ran
};
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};
HEJ::optional<HEJFOG::Unweighter> unweighter{};
std::map<HEJFOG::Status, int> status_counter;
std::vector<HEJ::Event> events;
int trials = 0;
// warm-up phase to train unweighter
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()];
if(generator.status() == HEJFOG::good && analysis->pass_cuts(ev, ev)) {
events.emplace_back(std::move(ev));
++warmup_progress;
}
}
std::cout << std::endl;
unweighter = HEJFOG::Unweighter{
begin(events), end(events), config.unweight->max_dev, *ran,
config.jets.peak_pt?(*config.jets.peak_pt):0.
};
std::vector<HEJ::Event> unweighted_events;
for(auto && ev: events) {
auto unweighted = unweighter->unweight(std::move(ev));
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";
}
HEJ::ProgressBar<long long> progress{std::cout, config.events};
progress.increment(events.size());
events.reserve(config.events);
for(; events.size() < static_cast<size_t>(config.events); ++trials){
auto ev = generator.gen_event();
++status_counter[generator.status()];
if(generator.status() == HEJFOG::good && analysis->pass_cuts(ev, ev)) {
if(unweighter) {
auto unweighted = unweighter->unweight(std::move(ev));
if(! unweighted) continue;
ev = std::move(*unweighted);
}
events.emplace_back(std::move(ev));
++progress;
}
}
std::cout << std::endl;
- HEJFOG::CrossSectionAccumulator xs;
+ HEJ::CrossSectionAccumulator xs;
for(auto & ev: events){
rescale_weights(ev, invGeV2_to_pb/trials);
analysis->fill(ev, ev);
writer.write(ev);
xs.fill(ev);
}
analysis->finalise();
const std::chrono::duration<double> run_time = (clock::now() - start_time);
std::cout << "\nTask Runtime: " << run_time.count() << " seconds.\n\n";
- std::cout.precision(10);
- std::cout.setf(std::ios::fixed);
- std::cout
- << " " << std::left << std::setw(20)
- << "Cross section: " << xs.total().value
- << " +- " << std::sqrt(xs.total().error) << " (pb)\n";
- for(auto const & xs_type: xs) {
- std::cout
- << " " << std::left << std::setw(20)
- << (HEJ::event_type::names[xs_type.first] + ": "s)
- << xs_type.second.value << " +- "
- << std::sqrt(xs_type.second.error) << " (pb)\n";
- }
+ std::cout << xs << '\n';
- std::cout << '\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 << "%\n";
}
}
diff --git a/include/HEJ/CrossSectionAccumulator.hh b/include/HEJ/CrossSectionAccumulator.hh
new file mode 100644
index 0000000..4d5b91e
--- /dev/null
+++ b/include/HEJ/CrossSectionAccumulator.hh
@@ -0,0 +1,68 @@
+#pragma once
+
+#include <map>
+#include <ostream>
+#include <string>
+
+#include "HEJ/Event.hh"
+#include "HEJ/event_types.hh"
+
+namespace HEJ {
+ template<typename T>
+ struct XSWithError {
+ T value = T{};
+ T error = T{};
+ };
+
+ /**
+ * @brief Sum of Cross Section for different subproccess
+ */
+ class CrossSectionAccumulator {
+ public:
+ void fill(HEJ::Event const & ev) {
+ const double wt = ev.central().weight;
+ auto & entry = xs_[ev.type()];
+ entry.value += wt;
+ entry.error += wt*wt;
+ total_.value += wt;
+ total_.error += wt*wt;
+ }
+
+ auto begin() const {
+ return std::begin(xs_);
+ }
+
+ auto end() const {
+ return std::end(xs_);
+ }
+
+ //! total Cross Section and error
+ XSWithError<double> total() const {
+ return total_;
+ }
+
+
+ private:
+ std::map<HEJ::event_type::EventType, XSWithError<double>> xs_;
+ XSWithError<double> total_;
+ };
+
+ std::ostream& operator<<(std::ostream& os, const CrossSectionAccumulator& xs){
+ std::streamsize ss = std::cout.precision();
+ os << std::scientific << std::setprecision(10)
+ << " " << std::left << std::setw(20)
+ << "Cross section: " << xs.total().value
+ << " +- " << std::sqrt(xs.total().error) << " (pb)\n";
+ for(auto const & xs_type: xs) {
+ os
+ << " " << std::left << std::setw(20)
+ << (HEJ::event_type::names[xs_type.first] + std::string(": "))
+ << xs_type.second.value << " +- "
+ << std::sqrt(xs_type.second.error) << " (pb)\n";
+ }
+ os << std::defaultfloat;
+ std::cout.precision(ss);
+ return os;
+ }
+
+}
diff --git a/src/bin/HEJ.cc b/src/bin/HEJ.cc
index ca81ec7..f0c672a 100644
--- a/src/bin/HEJ.cc
+++ b/src/bin/HEJ.cc
@@ -1,185 +1,191 @@
/**
* \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie
* \date 2019
* \copyright GPLv2 or later
*/
#include <array>
#include <chrono>
#include <iostream>
#include <limits>
#include <memory>
#include <numeric>
#include "yaml-cpp/yaml.h"
#include "LHEF/LHEF.h"
#include "fastjet/ClusterSequence.hh"
#include "HEJ/CombinedEventWriter.hh"
#include "HEJ/config.hh"
+#include "HEJ/CrossSectionAccumulator.hh"
#include "HEJ/Event.hh"
#include "HEJ/EventReweighter.hh"
#include "HEJ/get_analysis.hh"
#include "HEJ/make_RNG.hh"
#include "HEJ/ProgressBar.hh"
#include "HEJ/stream.hh"
#include "HEJ/Version.hh"
#include "HEJ/YAMLreader.hh"
int event_number(std::string const & record){
size_t start = record.rfind("Number of Events");
start = record.find_first_of("123456789", start);
if(start == std::string::npos) {
throw std::invalid_argument("no event number record found");
}
const size_t end = record.find_first_not_of("0123456789", start);
return std::stoi(record.substr(start, end - start));
}
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
){
try{
return HEJ::get_analysis(parameters);
}
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;
}
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;
}
std::cout << HEJ::Version::logo() << "\n" << std::endl;
fastjet::ClusterSequence::print_banner();
// read configuration
const HEJ::Config config = load_config(argv[1]);
HEJ::istream in{argv[2]};
std::unique_ptr<HEJ::Analysis> analysis = get_analysis(
config.analysis_parameters
);
assert(analysis != nullptr);
LHEF::Reader reader{in};
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();
HEJ::CombinedEventWriter writer{config.output, std::move(heprup)};
double global_reweight = 1.;
int max_events = std::numeric_limits<int>::max();
if(argn > 3){
max_events = std::stoi(argv[3]);
const int input_events = event_number(reader.headerBlock);
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);
assert(ran != nullptr);
HEJ::EventReweighter hej{
reader.heprup,
std::move(scale_gen),
to_EventReweighterConfig(config),
*ran
};
int 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;
// Loop over the events in the inputfile
while(reader.readEvent()){
// reweight events so that the total cross section is conserved
reader.hepeup.setWeight(0, global_reweight * reader.hepeup.weight());
if(nevent == max_events) break;
++nevent;
// calculate HEJ weight
HEJ::Event FO_event{
HEJ::UnclusteredEvent{reader.hepeup},
config.fixed_order_jets.def, config.fixed_order_jets.min_pt,
};
auto resummed_events = hej.reweight(FO_event, config.trials);
++nevent_type[FO_event.type()];
if(resummed_events.empty()) ++nfailed_type[FO_event.type()];
for(auto const & 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);
+ xs.fill(ev);
}
}
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 << '\n';
for(size_t ev_type = first_type; ev_type <= last_type; ++ev_type){
std::cout << '\t' << names[ev_type] << ": " << nevent_type[ev_type]
<< ", failed to reconstruct " << nfailed_type[ev_type]
<< '\n';
}
+
+ std::cout << '\n' << xs << '\n';
+
std::chrono::duration<double> run_time = (clock::now() - start_time);
std::cout << "Finished " << 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";
}

File Metadata

Mime Type
text/x-diff
Expires
Mon, Jan 20, 10:59 PM (1 d, 20 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4224613
Default Alt Text
(18 KB)

Event Timeline