Page MenuHomeHEPForge
Files F8723670

No OneTemporary
Actions

View Options

diff --git a/FixedOrderGen/include/EventFilter.hh b/FixedOrderGen/include/EventFilter.hh
deleted file mode 100644
index 5d5be65..0000000
--- a/FixedOrderGen/include/EventFilter.hh
+++ /dev/null
@@ -1,60 +0,0 @@
-#pragma once
-
-#include <vector>
-
-#include "RHEJ/Event.hh"
-
-namespace CLHEP {
- class Ranlux64Engine;
-}
-
-namespace HEJFOG{
- struct EventFilter {
- public:
- virtual EventFilter & operator>>(RHEJ::Event & ev) = 0;
- virtual EventFilter & operator<<(RHEJ::Event const & ev) = 0;
- virtual EventFilter & operator<<(RHEJ::Event && ev) = 0;
- virtual explicit operator bool() = 0;
- virtual ~EventFilter() = default;
- };
-
- class TrivialFilter: public EventFilter {
- public:
- TrivialFilter() = default;
- EventFilter & operator>>(RHEJ::Event & ev) override;
- EventFilter & operator<<(RHEJ::Event const & ev) override;
- EventFilter & operator<<(RHEJ::Event && ev) override;
- explicit operator bool() override;
- ~TrivialFilter() = default;
-
- private:
- std::vector<RHEJ::Event> store_;
- bool good_;
- };
-
- class Unweighter : public EventFilter {
- public:
- Unweighter(
- double max_dev, int sample_size,
- CLHEP::Ranlux64Engine & ran
- );
-
- EventFilter & operator>>(RHEJ::Event & ev) override;
- EventFilter & operator<<(RHEJ::Event const & ev) override;
- EventFilter & operator<<(RHEJ::Event && ev) override;
- explicit operator bool() override;
- ~Unweighter() override = default;
-
- private:
- void update_sample();
-
- std::vector<RHEJ::Event> sample_;
- bool sampling_finished_;
- bool good_;
- size_t sample_size_;
- double max_dev_;
- double cut_;
- CLHEP::Ranlux64Engine & ran_;
- };
-
-}
diff --git a/FixedOrderGen/src/EventFilter.cc b/FixedOrderGen/src/EventFilter.cc
deleted file mode 100644
index 86b4fa8..0000000
--- a/FixedOrderGen/src/EventFilter.cc
+++ /dev/null
@@ -1,120 +0,0 @@
-#include "EventFilter.hh"
-
-#include <cassert>
-
-#include <CLHEP/Random/Randomize.h>
-#include <CLHEP/Random/RanluxEngine.h>
-
-namespace HEJFOG{
-
- EventFilter & TrivialFilter::operator<<(RHEJ::Event const & ev) {
- store_.emplace_back(ev);
- good_ = true;
- return *this;
- }
-
- EventFilter & TrivialFilter::operator<<(RHEJ::Event && ev) {
- store_.emplace_back(std::move(ev));
- good_ = true;
- return *this;
- }
-
- EventFilter & TrivialFilter::operator>>(RHEJ::Event & ev){
- if(! store_.empty()) {
- ev = std::move(store_.back());
- store_.pop_back();
- }
- else good_ = false;
- return *this;
- }
-
- TrivialFilter::operator bool(){
- return good_;
- }
-
- Unweighter::Unweighter(
- double max_dev, int sample_size,
- CLHEP::Ranlux64Engine & ran
- ):
- sampling_finished_{false},
- good_{false},
- sample_size_{static_cast<size_t>(sample_size)},
- max_dev_{max_dev},
- ran_{ran}
- {
- assert(sample_size > 0);
- sample_.reserve(sample_size);
- }
-
- EventFilter & Unweighter::operator<<(RHEJ::Event const & ev) {
- sample_.emplace_back(ev);
- update_sample();
- good_ = sampling_finished_;
- return *this;
- }
-
- EventFilter & Unweighter::operator<<(RHEJ::Event && ev){
- sample_.emplace_back(std::move(ev));
- update_sample();
- good_ = sampling_finished_;
- return *this;
- }
-
- namespace {
- template<typename Iterator>
- std::pair<double, double> logweight_mean_err(Iterator begin, Iterator end){
- double mean = 0.;
- double err = 0.;
- double awt_sum = 0.;
- do{
- const double awt = std::abs(begin->central().weight);
- const double tmp = awt*std::log(awt);
- mean += tmp;
- err += tmp*tmp;
- awt_sum += awt;
- } while(++begin != end);
- mean /= awt_sum;
- err = std::sqrt(err)/awt_sum;
- return std::make_pair(mean, err);
- }
- }
-
- void Unweighter::update_sample(){
- if(sampling_finished_) return;
- if(sample_.size() >= sample_size_){
- sampling_finished_ = true;
- double mean, err;
- std::tie(mean, err) = logweight_mean_err(begin(sample_), end(sample_));
- cut_ = std::exp(mean + max_dev_*err);
- }
- }
-
- namespace {
- void normalise_weights(RHEJ::Event & ev, double target) {
- const double awt = std::abs(ev.central().weight);
- ev.central().weight *= target/awt;
- for(auto & var: ev.variations()) var.weight *= target/awt;
- }
- }
-
- EventFilter & Unweighter::operator>>(RHEJ::Event & ev){
- if(!sampling_finished_) return *this;
- while(!sample_.empty()) {
- auto cur_ev = std::move(sample_.back());
- sample_.pop_back();
- const double awt = std::abs(cur_ev.central().weight);
- if(ran_.flat() < awt/cut_) {
- if(awt < cut_) normalise_weights(cur_ev, cut_);
- ev = std::move(cur_ev);
- return *this;
- }
- }
- good_ = false;
- return *this;
- }
-
- Unweighter::operator bool(){
- return good_;
- }
-
-}
diff --git a/FixedOrderGen/src/main.cc b/FixedOrderGen/src/main.cc
index 93fc71d..4f64f84 100644
--- a/FixedOrderGen/src/main.cc
+++ b/FixedOrderGen/src/main.cc
@@ -1,205 +1,192 @@
/**
* Name: main.cc
* Authors: Jeppe R. Andersen
*/
#include <fstream>
#include <algorithm>
#include <memory>
#include <chrono>
#include <iostream>
#include <map>
#include "yaml-cpp/yaml.h"
#include "config.hh"
#include "LHEF/LHEF.h"
#include "RHEJ/CombinedEventWriter.hh"
#include "RHEJ/get_analysis.hh"
#include "RHEJ/utility.hh"
//#include "RHEJ/EventReweighter.hh"
#include "RHEJ/stream.hh"
#include "RHEJ/LesHouchesWriter.hh"
#include "EventGenerator.hh"
#include "PhaseSpacePoint.hh"
#include "Ranlux64Engine.hh"
-#include "EventFilter.hh"
namespace{
constexpr auto banner =
" __ ___ __ ______ __ __ \n"
" / / / (_)___ _/ /_ / ____/___ ___ _________ ___ __ / /__ / /______ \n"
" / /_/ / / __ `/ __ \\ / __/ / __ \\/ _ \\/ ___/ __ `/ / / / __ / / _ \\/ __/ ___/ \n"
" / __ / / /_/ / / / / / /___/ / / / __/ / / /_/ / /_/ / / /_/ / __/ /_(__ ) \n"
" /_/ /_/_/\\__, /_/ /_/ /_____/_/ /_/\\___/_/ \\__, /\\__, / \\____/\\___/\\__/____/ \n"
" ____///__/ __ ____ ///__//____/ ______ __ \n"
" / ____(_) _____ ____/ / / __ \\_________/ /__ _____ / ____/__ ____ ___ _________ _/ /_____ _____\n"
" / /_ / / |/_/ _ \\/ __ / / / / / ___/ __ / _ \\/ ___/ / / __/ _ \\/ __ \\/ _ \\/ ___/ __ `/ __/ __ \\/ ___/\n"
" / __/ / /> </ __/ /_/ / / /_/ / / / /_/ / __/ / / /_/ / __/ / / / __/ / / /_/ / /_/ /_/ / / \n"
" /_/ /_/_/|_|\\___/\\__,_/ \\____/_/ \\__,_/\\___/_/ \\____/\\___/_/ /_/\\___/_/ \\__,_/\\__/\\____/_/ \n"
;
constexpr double invGeV2_to_pb = 389379292.;
}
HEJFOG::Config load_config(char const * filename){
try{
return HEJFOG::load_config(filename);
}
catch(std::exception const & exc){
std::cerr << "Error: " << exc.what() << '\n';
std::exit(EXIT_FAILURE);
}
}
std::unique_ptr<RHEJ::Analysis> get_analysis(
YAML::Node const & parameters
){
try{
return RHEJ::get_analysis(parameters);
}
catch(std::exception const & exc){
std::cerr << "Failed to load analysis: " << exc.what() << '\n';
std::exit(EXIT_FAILURE);
}
}
std::string progress_bar(){
std::string result = "0% ";
for(int i = 10; i <= 100; i+= 10){
result += " " + std::to_string(i) + "%";
}
result += "\n|";
for(int i = 10; i <= 100; i+= 10){
result += "---------|";
}
return result + '\n';
}
int main(int argn, char** argv) {
if (argn < 2) {
std::cerr << "\n# Usage:\n.FOgen config_file\n";
return EXIT_FAILURE;
}
std::cout << banner;
fastjet::ClusterSequence::print_banner();
using clock = std::chrono::system_clock;
const auto start_time = clock::now();
// read configuration
auto config = load_config(argv[1]);
std::unique_ptr<RHEJ::Analysis> analysis = get_analysis(
config.analysis_parameters
);
assert(analysis != nullptr);
auto ran = HEJFOG::make_Ranlux64Engine();
if(config.RanLux_init) ran.restoreStatus(config.RanLux_init->c_str());
HEJFOG::EventGenerator generator{
config.process,
config.beam,
config.scale_gen,
config.jets,
config.pdf_id,
config.unordered_fraction,
config.Higgs_properties,
config.Higgs_coupling,
ran
};
//TODO: fix Les Houches init block (HEPRUP)
LHEF::HEPRUP lhefheprup;
lhefheprup.IDBMUP=std::pair<long,long>(config.beam.particles[0], config.beam.particles[1]);
lhefheprup.EBMUP=std::make_pair(config.beam.energy, config.beam.energy);
lhefheprup.PDFGUP=std::make_pair(0,0);
lhefheprup.PDFSUP=std::make_pair(config.pdf_id,config.pdf_id);
lhefheprup.NPRUP=1;
lhefheprup.XSECUP=std::vector<double>(1.);
lhefheprup.XERRUP=std::vector<double>(1.);
lhefheprup.LPRUP=std::vector<int>{1};
RHEJ::CombinedEventWriter writer{config.output, lhefheprup};
std::map<HEJFOG::Status, int> status_counter;
- std::unique_ptr<HEJFOG::EventFilter> filter =
- RHEJ::make_unique<HEJFOG::TrivialFilter>();
- if(config.unweight) {
- filter = RHEJ::make_unique<HEJFOG::Unweighter>(
- config.unweight->max_dev, config.unweight->sample_size, ran
- );
- }
-
// Perform N trial generations
int iprint = 0;
double FKL_xs = 0., FKL_xs_err = 0.;
double uno_xs = 0., uno_xs_err = 0.;
std::cout << '\n' << progress_bar();
for (int trials = 0; trials< config.trials;trials++){
if (trials==iprint) {
std::cout << ".";
std::cout.flush();
iprint+=(int) config.trials/100;
}
auto ev = generator.gen_event();
++status_counter[generator.status()];
if(generator.status() != HEJFOG::good) continue;
ev.central().weight /= config.trials;
ev.central().weight *= invGeV2_to_pb;
for(auto & var: ev.variations()){
var.weight /= config.trials;
var.weight *= invGeV2_to_pb;
}
+ if(analysis->pass_cuts(ev, ev)){
+ analysis->fill(ev, ev);
+ writer.write(ev);
- *filter << std::move(ev);
- while(*filter >> ev) {
- if(analysis->pass_cuts(ev, ev)){
- analysis->fill(ev, ev);
- writer.write(ev);
-
- if(ev.type() == RHEJ::event_type::FKL){
- FKL_xs += ev.central().weight;
- FKL_xs_err += ev.central().weight*ev.central().weight;
- }
- else {
- uno_xs += ev.central().weight;
- uno_xs_err += ev.central().weight*ev.central().weight;
- }
+ if(ev.type() == RHEJ::event_type::FKL){
+ FKL_xs += ev.central().weight;
+ FKL_xs_err += ev.central().weight*ev.central().weight;
+ }
+ else {
+ uno_xs += ev.central().weight;
+ uno_xs_err += ev.central().weight*ev.central().weight;
}
}
} // main event loop
std::cout << std::endl;
std::chrono::duration<double> run_time = (clock::now() - start_time);
std::cout << "\nTask Runtime: " << run_time.count() << " seconds.\n";
std::cout.precision(10);
std::cout.width(12);
std::cout.setf(std::ios::fixed);
std::cout
<< "\nCross section: " << (FKL_xs + uno_xs)
<< " +- " << std::sqrt(FKL_xs_err + uno_xs_err) << " (pb)"
"\n FKL: " << FKL_xs << " +- " << std::sqrt(FKL_xs_err) << " (pb)"
"\n unordered: " << uno_xs << " +- " << std::sqrt(uno_xs_err) << " (pb)"
"\n";
std::cout << '\n';
for(auto && entry: status_counter){
const double fraction = (double) entry.second/config.trials;
const int percent = std::round(100*fraction);
std::cout << "status "
<< std::left << std::setw(16) << (to_string(entry.first) + ":")
<< " [";
for(int i = 0; i < percent/2; ++i) std::cout << '#';
for(int i = percent/2; i < 50; ++i) std::cout << ' ';
std::cout << "] " << percent << "%\n";
}
}

File Metadata

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

Event Timeline