Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F8723697
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
5 KB
Subscribers
None
View Options
diff --git a/src/main.cc b/src/main.cc
index 19c1fba..f86693d 100644
--- a/src/main.cc
+++ b/src/main.cc
@@ -1,163 +1,176 @@
/**
* Name: main.cc
* Authors: Tuomas Hapola, Andreas Maier <andreas.maier@durham.ac.uk>
*
*/
#include <fstream>
#include <algorithm>
#include <memory>
#include <chrono>
#include <iostream>
#include "yaml-cpp/yaml.h"
#include "RHEJ/CombinedEventWriter.hh"
#include "RHEJ/config.hh"
#include "RHEJ/EventReweighter.hh"
#include "RHEJ/get_analysis.hh"
#include "LHEF/LHEF.h"
#include "RHEJ/utility.hh"
#include "RHEJ/Version.hh"
#include "RHEJ/stream.hh"
#include "RHEJ/YAMLreader.hh"
+#include "RHEJ/ProgressBar.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));
}
RHEJ::Config load_config(char const * filename){
try{
return RHEJ::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 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 " << RHEJ::Version::package_name_full()
<< ", revision " << RHEJ::Version::revision() << " ("
<< time_to_string(clock::to_time_t(start_time)) << ")" << std::endl;
}
+ fastjet::ClusterSequence::print_banner();
// read configuration
const RHEJ::Config config = load_config(argv[1]);
RHEJ::istream in{argv[2]};
std::unique_ptr<RHEJ::Analysis> analysis = get_analysis(
config.analysis_parameters
);
assert(analysis != nullptr);
LHEF::Reader reader{in};
RHEJ::CombinedEventWriter writer{config.output, reader.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";
}
RHEJ::ScaleGenerator scale_gen{
config.scales.base,
config.scales.factors,
config.scales.max_ratio
};
RHEJ::EventReweighter rhej{
reader.heprup,
std::move(scale_gen),
to_EventReweighterConfig(config)
};
if(config.RanLux_init){
RHEJ::PhaseSpacePoint::reset_ranlux(*config.RanLux_init);
}
int nevent = 0;
std::array<int, RHEJ::event_type::last_type + 1>
nevent_type{0}, nfailed_type{0};
+ // unique_ptr is a workaround:
+ // RHEJ::optional is a better fit, but gives spurious errors with g++ 7.3.0
+ std::unique_ptr<RHEJ::ProgressBar<double>> progress = nullptr;
+ if(! reader.heprup.XSECUP.empty()) {
+ const double Born_xs = std::accumulate(
+ begin(reader.heprup.XSECUP), end(reader.heprup.XSECUP), 0.
+ );
+ progress = std::make_unique<RHEJ::ProgressBar<double>>(std::cout, Born_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;
if (nevent % 10000 == 0){
std::cout << "Passed " << nevent << " events ("
<< time_to_string(clock::to_time_t(clock::now())) << ")"<< std::endl;
}
// calculate rHEJ weight
RHEJ::Event FO_event{
RHEJ::UnclusteredEvent{reader.hepeup},
config.fixed_order_jets.def, config.fixed_order_jets.min_pt,
};
auto resummed_events = rhej.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);
}
}
+ if(progress) progress->increment(FO_event.central().weight);
} // main event loop
+ std::cout << '\n';
using namespace RHEJ::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::chrono::duration<double> run_time = (clock::now() - start_time);
std::cout << "Finished " << RHEJ::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
Details
Attached
Mime Type
text/x-diff
Expires
Mon, Jan 20, 9:12 PM (22 h, 3 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4242409
Default Alt Text
(5 KB)
Attached To
rHEJ HEJ
Event Timeline
Log In to Comment