diff --git a/t/test_ME_generic.cc b/t/test_ME_generic.cc index 5d29742..16ebc54 100644 --- a/t/test_ME_generic.cc +++ b/t/test_ME_generic.cc @@ -1,141 +1,147 @@ /** * \brief Generic tester for the ME for a given set of PSP * * \note reference weights and PSP (as LHE file) have to be given as * _individual_ files * * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019 * \copyright GPLv2 or later */ +#include <algorithm> #include <fstream> +#include <math.h> #include <random> -#include <algorithm> #include "LHEF/LHEF.h" -#include "HEJ/MatrixElement.hh" #include "HEJ/Event.hh" -#include "HEJ/YAMLreader.hh" +#include "HEJ/MatrixElement.hh" #include "HEJ/stream.hh" +#include "HEJ/YAMLreader.hh" constexpr double alpha_s = 0.118; constexpr double ep = 1e-5; void shuffle_particles(HEJ::Event::EventData & ev) { static std::mt19937_64 ran{0}; std::shuffle(begin(ev.incoming), end(ev.incoming), ran); std::shuffle(begin(ev.outgoing), end(ev.outgoing), ran); } void dump(HEJ::Event const & ev){ { LHEF::Writer writer{std::cout}; std::cout << std::setprecision(6); writer.hepeup = to_HEPEUP(std::move(ev), nullptr); writer.writeEvent(); } std::cout << "Rapidity ordering:\n"; for(const auto & part: ev.outgoing()){ std::cout << std::setw(2) << part.type << ": "<< std::setw(7) << part.rapidity() << std::endl; } } enum MEComponent {tree, virt}; MEComponent guess_component(std::string const & data_file) { if(data_file.find("virt") != data_file.npos) return MEComponent::virt; return MEComponent::tree; } int main(int argn, char** argv){ if(argn != 4 && argn != 5){ std::cerr << "\n# Usage:\n."<< argv[0] <<" config.yml ME_weights input_file.lhe\n\n"; return EXIT_FAILURE; } bool OUTPUT_MODE = false; if(argn == 5 && std::string("OUTPUT")==std::string(argv[4])) OUTPUT_MODE = true; const HEJ::Config config = HEJ::load_config(argv[1]); std::fstream wgt_file; if ( OUTPUT_MODE ) { std::cout << "_______________________USING OUTPUT MODE!_______________________" << std::endl; wgt_file.open(argv[2], std::fstream::out); wgt_file.precision(10); } else { wgt_file.open(argv[2], std::fstream::in); } HEJ::istream in{argv[3]}; LHEF::Reader reader{in}; const auto component = guess_component(argv[2]); HEJ::MatrixElement ME{ [](double){ return alpha_s; }, HEJ::to_MatrixElementConfig(config) }; double max_ratio = 0.; size_t idx_max_ratio = 0; HEJ::Event ev_max_ratio(HEJ::Event::EventData{}.cluster( config.resummation_jets.def,0 ) ); double av_ratio = 0; size_t i = 0; while(reader.readEvent()){ ++i; HEJ::Event::EventData data{reader.hepeup}; shuffle_particles(data); HEJ::Event event{ data.cluster( config.resummation_jets.def, config.resummation_jets.min_pt ) }; const double our_ME = (component == MEComponent::tree)? ME.tree(event).central: ME.virtual_corrections(event).central ; + if(!isfinite(our_ME)){ + std::cerr << "Found non-finite ME ("<< our_ME <<")\n"; + dump(event); + return EXIT_FAILURE; + } if ( OUTPUT_MODE ) { wgt_file << our_ME << std::endl; } else { std::string line; if(!std::getline(wgt_file,line)) break; const double ref_ME = std::stod(line); const double diff = std::abs(our_ME/ref_ME-1.); av_ratio+=diff; if( diff > max_ratio ) { max_ratio = diff; idx_max_ratio = i; ev_max_ratio = event; } if( diff > ep ){ size_t precision(std::cout.precision()); std::cout.precision(16); std::cout<< "Large difference in PSP " << i << "\nis: "<<our_ME << " should: " << ref_ME << " => difference: " << diff << std::endl; std::cout.precision(precision); dump(event); return EXIT_FAILURE; } } } wgt_file.close(); if ( i<100 ) throw std::invalid_argument{"Not enough PSP tested"}; if ( !OUTPUT_MODE ) { size_t precision(std::cout.precision()); std::cout.precision(16); std::cout << "Avg ratio after " << i << " PSP: " << av_ratio/i << std::endl; std::cout << "maximal ratio at " << idx_max_ratio << ": " << max_ratio << std::endl; std::cout.precision(precision); } return EXIT_SUCCESS; }