diff --git a/src/HepMCWriter.cc b/src/HepMCWriter.cc index 3fe2817..4aca4b3 100644 --- a/src/HepMCWriter.cc +++ b/src/HepMCWriter.cc @@ -1,243 +1,243 @@ #include "RHEJ/HepMCWriter.hh" #include <cassert> #ifdef RHEJ_BUILD_WITH_HepMC_VERSION #if RHEJ_BUILD_WITH_HepMC_VERSION >= 3 #include "HepMC/WriterAscii.h" #include "HepMC/LHEFAttributes.h" #else #include "HepMC/IO_GenEvent.h" #endif #include "HepMC/GenVertex.h" #include "HepMC/GenParticle.h" #include "RHEJ/config.hh" #include "RHEJ/utility.hh" #include "RHEJ/Event.hh" #include "RHEJ/exceptions.hh" #include "RHEJ/Version.hh" namespace RHEJ{ #if RHEJ_BUILD_WITH_HepMC_VERSION >= 3 namespace { void add_generator_tag(LHEF::HEPRUP & heprup){ heprup.generators.emplace_back(LHEF::XMLTag{}); - heprup.generators.back().name = Version().package_name(); - heprup.generators.back().version = Version().String(); + heprup.generators.back().name = Version::package_name(); + heprup.generators.back().version = Version::String(); } void reset_weight_info(LHEF::HEPRUP & heprup){ heprup.IDWTUP = -4; // use placeholders for unknown init block values // we can overwrite them after processing all events heprup.XSECUP = {0.}; heprup.XERRUP = {0.}; heprup.XMAXUP = {0.}; } HepMC::FourVector to_FourVector(Sparticle const & sp){ return {sp.px(), sp.py(), sp.pz(), sp.E()}; } HepMC::shared_ptr<HepMC::GenRunInfo> & add_HEPRUP( HepMC::shared_ptr<HepMC::GenRunInfo> & runinfo, LHEF::HEPRUP heprup ){ add_generator_tag(heprup); reset_weight_info(heprup); auto hepr = HepMC::make_shared<HepMC::HEPRUPAttribute>(); hepr->heprup = std::move(heprup); runinfo->add_attribute("HEPRUP", hepr); for (int i = 0, N = hepr->heprup.generators.size(); i < N; ++i ){ HepMC::GenRunInfo::ToolInfo tool; tool.name = hepr->heprup.generators[i].name; tool.version = hepr->heprup.generators[i].version; tool.description = hepr->heprup.generators[i].contents; runinfo->tools().push_back(tool); } return runinfo; } } // namespace anonymous struct HepMCWriter::HepMCWriterImpl: public EventWriter{ Config config_; HepMC::shared_ptr<HepMC::GenRunInfo> runinfo_; HepMC::WriterAscii writer_; HepMCWriterImpl( std::string const & file, Config const & conf, LHEF::HEPRUP && heprup ): config_(conf), runinfo_{HepMC::make_shared<HepMC::GenRunInfo>()}, writer_{file, add_HEPRUP(runinfo_, std::move(heprup))} {} HepMCWriterImpl & operator=(HepMCWriterImpl const & other) = delete; HepMCWriterImpl(HepMCWriterImpl const & other) = delete; HepMCWriterImpl & operator=(HepMCWriterImpl && other) = delete; HepMCWriterImpl(HepMCWriterImpl && other) = delete; ~HepMCWriterImpl(){ writer_.close(); } void write(Event const & ev){ if(!ev.decays().empty()){ throw not_implemented{"HepMC output for events with decays"}; } auto vx = HepMC::make_shared<HepMC::GenVertex>(); for(auto const & in: ev.incoming()){ vx->add_particle_in( HepMC::make_shared<HepMC::GenParticle>( to_FourVector(in), static_cast<int>(in.type), -1 ) ); } for(auto const & out: ev.outgoing()){ vx->add_particle_out( HepMC::make_shared<HepMC::GenParticle>( to_FourVector(out), static_cast<int>(out.type), +1 ) ); } HepMC::GenEvent out_ev{runinfo_, HepMC::Units::GEV, HepMC::Units::MM}; out_ev.add_vertex(vx); out_ev.weights().reserve(ev.variations().size() + 1u); out_ev.weights().emplace_back(ev.central().weight); for(auto const & var: ev.variations()){ out_ev.weights().emplace_back(var.weight); } writer_.write_event(out_ev); } }; void HepMCWriter::HepMCWriterImplDeleter::operator()(HepMCWriterImpl* p) { delete p; } #else // HepMC 2 namespace{ HepMC::FourVector to_FourVector(Sparticle const & sp){ return {sp.px(), sp.py(), sp.pz(), sp.E()}; } } struct HepMCWriter::HepMCWriterImpl: public EventWriter{ struct HepMCInfo{ size_t event_count = 0.; double tot_weight = 0.; double tot_weight2 = 0.; void add_event(Event const & ev){ double wt = ev.central().weight; tot_weight += wt; tot_weight2 += wt * wt; ++event_count; } }; Config config_; HepMC::IO_GenEvent writer_; HepMCInfo genInfo_; HepMCWriterImpl( std::string const & file, Config const & conf, LHEF::HEPRUP && ): config_(conf), writer_{file} {} HepMCWriterImpl & operator=(HepMCWriterImpl const & other) = delete; HepMCWriterImpl(HepMCWriterImpl const & other) = delete; HepMCWriterImpl & operator=(HepMCWriterImpl && other) = delete; HepMCWriterImpl(HepMCWriterImpl && other) = delete; void write(Event const & ev){ if(!ev.decays().empty()){ throw not_implemented{"HepMC output for events with decays"}; } auto vx = new HepMC::GenVertex(); for(auto const & in: ev.incoming()){ vx->add_particle_in( new HepMC::GenParticle{ to_FourVector(in), static_cast<int>(in.type), -1 } ); } for(auto const & out: ev.outgoing()){ vx->add_particle_out( new HepMC::GenParticle{ to_FourVector(out), static_cast<int>(out.type), +1 } ); } HepMC::GenEvent out_ev{HepMC::Units::GEV, HepMC::Units::MM}; out_ev.add_vertex(vx); /// weights out_ev.weights().push_back(ev.central().weight); for(auto const & var: ev.variations()){ out_ev.weights().push_back(var.weight); } /// @TODO add name list /// general informations genInfo_.add_event(ev); out_ev.set_signal_process_id(ev.type()+1); /// "+1": conistent with lhe out_ev.set_event_scale(ev.central().mur); /// event scale out_ev.set_event_number(genInfo_.event_count); /// event number /// cross section HepMC::GenCrossSection t_xs; t_xs.set_cross_section( genInfo_.tot_weight , sqrt(genInfo_.tot_weight2) ); out_ev.set_cross_section( t_xs ); /// @TODO add alphaQCD (need function) and alphaQED /// @TODO output pdf (currently not avaiable from event alone) writer_.write_event(&out_ev); } }; void HepMCWriter::HepMCWriterImplDeleter::operator()(HepMCWriterImpl* p) { delete p; } #endif // HepMC 2 HepMCWriter::HepMCWriter(std::string const & file, Config const & conf, LHEF::HEPRUP heprup): impl_{std::unique_ptr<HepMCWriterImpl, HepMCWriterImplDeleter>{ new HepMCWriterImpl(file, conf, std::move(heprup)) }} {} void HepMCWriter::write(Event const & ev){ impl_->write(ev); } } // namespace RHEJ #else // no HepMC namespace RHEJ{ class HepMCWriter::HepMCWriterImpl{}; HepMCWriter::HepMCWriter(std::string const &, Config const & conf, LHEF::HEPRUP){ throw std::invalid_argument( "Failed to create HepMC writer: " "Reversed HEJ was built without HepMC support" ); } void HepMCWriter::write(Event const &){ assert(false); } void HepMCWriter::HepMCWriterImplDeleter::operator()(HepMCWriterImpl* p) { delete p; } } #endif diff --git a/src/LesHouchesWriter.cc b/src/LesHouchesWriter.cc index 57a4726..7561dad 100644 --- a/src/LesHouchesWriter.cc +++ b/src/LesHouchesWriter.cc @@ -1,92 +1,92 @@ #include <stdexcept> #include <memory> #include <cassert> #include "RHEJ/LesHouchesWriter.hh" #include "RHEJ/event_types.hh" #include "RHEJ/config.hh" #include "RHEJ/Event.hh" #include "RHEJ/Version.hh" namespace RHEJ{ namespace{ template<class T, class... Args> std::unique_ptr<T> make_unique(Args&&... a){ return std::unique_ptr<T>{new T{std::forward<Args>(a)...}}; } } LesHouchesWriter::LesHouchesWriter( std::string const & file, Config const & conf, LHEF::HEPRUP heprup ): config_(conf), out_{file, std::fstream::in | std::fstream::out | std::fstream::trunc}, writer_{RHEJ::make_unique<LHEF::Writer>(out_)} { if(! out_.is_open()){ throw std::ios_base::failure("Failed to open " + file); }; writer_->heprup = std::move(heprup); // lhe Stardard: IDWTUP (negative => weights = +/-) // 3: weight=+/-, xs given in head (same as default MG) // 4: weight=+/-, xs = avg(weights) writer_->heprup.IDWTUP = -3; writer_->heprup.generators.emplace_back(LHEF::XMLTag{}); - writer_->heprup.generators.back().name = Version().package_name(); - writer_->heprup.generators.back().version = Version().String(); + writer_->heprup.generators.back().name = Version::package_name(); + writer_->heprup.generators.back().version = Version::String(); // use placeholders for unknown init block values // we can overwrite them after processing all events writer_->heprup.XSECUP = std::vector<double>(event_type::last_type+1, 0.); writer_->heprup.XERRUP = std::vector<double>(event_type::last_type+1, 0.); writer_->heprup.XMAXUP = std::vector<double>(event_type::last_type+1, 0.); write_init(); } void LesHouchesWriter::write(Event const & ev){ assert(writer_ && out_.is_open()); const double wt = ev.central().weight; writer_->hepeup = RHEJ::to_HEPEUP(std::move(ev), &heprup()); writer_->writeEvent(); heprup().XSECUP[ev.type()] += wt; heprup().XERRUP[ev.type()] += wt*wt; if(wt > heprup().XMAXUP[ev.type()]){ heprup().XMAXUP[ev.type()] = wt; } } // this function is called after overwritting the Les Houches init block // assert that we have overwritten *exactly* the init block, // i.e. an intact event block is next void assert_next_event_intact(std::istream & out){ (void) out; // suppress compiler warnings if not in debug mode #ifndef NDEBUG std::string line; getline(out, line); assert(line == "<event>"); #endif } void LesHouchesWriter::rewrite_init(){ assert(writer_ && out_.is_open()); // replace placeholder entries const auto pos = out_.tellp(); out_.seekp(0); writer_->init(); assert_next_event_intact(out_); out_.seekp(pos); } LesHouchesWriter::~LesHouchesWriter(){ assert(writer_ && out_.is_open()); for(auto & xs_err: heprup().XERRUP) { xs_err = sqrt(xs_err); } rewrite_init(); } } diff --git a/src/main.cc b/src/main.cc index fd90541..4f034da 100644 --- a/src/main.cc +++ b/src/main.cc @@ -1,160 +1,160 @@ /** * 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" 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() << " (" + std::cout << "Starting " << RHEJ::Version::package_name_full() + << ", revision " << RHEJ::Version::revision() << " (" << time_to_string(clock::to_time_t(start_time)) << ")" << std::endl; } // 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, 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::EventReweighter rhej{ reader.heprup, 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}; // 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, config.scale_gen ); ++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); } } } // main event loop 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 " + 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"; }