Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F8723675
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
8 KB
Subscribers
None
View Options
diff --git a/FixedOrderGen/src/main.cc b/FixedOrderGen/src/main.cc
index c631673..46e0a46 100644
--- a/FixedOrderGen/src/main.cc
+++ b/FixedOrderGen/src/main.cc
@@ -1,215 +1,222 @@
/**
* 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 "ProgressBar.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);
}
}
void rescale_weights(RHEJ::Event & ev, double factor) {
ev.central().weight *= factor;
for(auto & var: ev.variations()){
var.weight *= factor;
}
}
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;
const int test_trials = config.trials;
double FKL_xs = 0., FKL_xs_err = 0.;
double uno_xs = 0., uno_xs_err = 0.;
std::vector<RHEJ::Event> events;
// warm-up phase: check how efficient we are
std::cout << "Warm-up ...\n";
HEJFOG::ProgressBar<int> warmup_progress{std::cout, test_trials};
const auto warmup_start = clock::now();
for(int trial = 0; trial < test_trials; ++trial){
++warmup_progress;
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));
}
}
std::cout << std::endl;
const double efficiency = static_cast<double>(events.size())/test_trials;
assert(0. <= efficiency && efficiency <= 1.);
std::cout << "Efficiency: " << efficiency << '\n';
- const int total_trials = test_trials/efficiency;
+ constexpr auto max_trials = std::numeric_limits<long long>::max();
+ if(max_trials*efficiency < test_trials) {
+ std::cerr << "Number of required configurations exceeds " << max_trials
+ << "\nPlease reduce the number of requested events.\n";
+ std::cout.setstate(std::ios_base::badbit); // prevent LHAPDF output
+ return EXIT_FAILURE;
+ }
+ const long long total_trials = test_trials/efficiency;
for(auto & ev: events) {
rescale_weights(ev, invGeV2_to_pb/total_trials);
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;
}
}
const auto warmup_end = clock::now();
const std::chrono::duration<double> warmup_time = warmup_end - warmup_start;
std::cout
<< "Generating " << (total_trials - test_trials) << " additional configurations\n"
"Estimated remaining time: "
<< warmup_time.count()/efficiency << " seconds.\n\n";
- HEJFOG::ProgressBar<int> progress{std::cout, total_trials};
+ HEJFOG::ProgressBar<long long> progress{std::cout, total_trials};
progress.increment(test_trials);
for(int trial = test_trials - 1; trial < total_trials; ++trial){
++progress;
auto ev = generator.gen_event();
++status_counter[generator.status()];
if(generator.status() == HEJFOG::good && analysis->pass_cuts(ev, ev)) {
rescale_weights(ev, invGeV2_to_pb/total_trials);
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;
}
}
}
std::cout << std::endl;
const 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 = static_cast<double>(entry.second)/total_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
Details
Attached
Mime Type
text/x-diff
Expires
Mon, Jan 20, 9:10 PM (23 h, 18 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4242389
Default Alt Text
(8 KB)
Attached To
rHEJ HEJ
Event Timeline
Log In to Comment