test_ME_generic.cc
View Options

	/**
	* \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 <cmath>
	#include <fstream>

	#include "HEJ/Event.hh"
	#include "HEJ/EventReader.hh"
	#include "HEJ/MatrixElement.hh"
	#include "HEJ/stream.hh"
	#include "HEJ/YAMLreader.hh"

	#include "hej_test.hh"

	constexpr double alpha_s = 0.118;
	constexpr double ep = 1e-5;
	constexpr double ep_mirror = 1e-3;

	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;
	}

	HEJ::Event::EventData mirror_event(HEJ::Event::EventData ev){
	for(auto & part: ev.incoming){
	auto & p{ part.p };
	p.reset(p.px(),p.py(),-p.pz(),p.E());
	}
	for(auto & part: ev.outgoing){
	auto & p{ part.p };
	p.reset(p.px(),p.py(),-p.pz(),p.E());
	}
	for(auto & decay: ev.decays){
	for(auto & part: decay.second){
	auto & p{ part.p };
	p.reset(p.px(),p.py(),-p.pz(),p.E());
	}
	}
	return ev;
	}

	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);
	}

	auto reader{ HEJ::make_reader(argv[3])};

	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->read_event()){
	++i;

	HEJ::Event::EventData data{reader->hepeup()};
	shuffle_particles(data);

	HEJ::Event::EventData data_mirror{mirror_event(data)};
	shuffle_particles(data_mirror);

	HEJ::Event event{
	data.cluster(
	config.resummation_jets.def,
	config.resummation_jets.min_pt
	)
	};
	HEJ::Event event_mirror{
	data_mirror.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(!std::isfinite(our_ME)){
	std::cerr << "Found non-finite ME ("<< our_ME <<")\n" << event << std::endl;
	return EXIT_FAILURE;
	}
	const double ME_mirror = (component == MEComponent::tree)?
	ME.tree(event_mirror).central:
	ME.virtual_corrections(event_mirror).central
	;
	if(!std::isfinite(ME_mirror)){
	std::cerr << "Found non-finite ME ("<< ME_mirror <<")\n" << event_mirror << std::endl;
	return EXIT_FAILURE;
	}

	if(std::abs(our_ME/ME_mirror-1.)>ep_mirror){
	size_t precision(std::cout.precision());
	std::cerr.precision(16);
	std::cerr<< "z-Mirrored ME gives different result " << i << "\n"
	<<our_ME << " vs " << ME_mirror << " => difference: "
	<< std::abs(our_ME/ME_mirror-1.) << "\n" << event
	<< "\nmirrored:\n" << event_mirror << std::endl;
	std::cerr.precision(precision);
	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::cerr.precision(16);
	std::cerr<< "Large difference in PSP " << i << "\nis: "<<our_ME
	<< " should: " << ref_ME << " => difference: " << diff << "\n"
	<< event << std::endl;
	std::cerr.precision(precision);
	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;
	}

File Metadata

Mime Type: text/x-c
Expires: Tue, Sep 30, 6:01 AM (1 d, 14 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 6533727
Default Alt Text: test_ME_generic.cc (4 KB)

test_ME_generic.cc
No OneTemporary
Actions

test_ME_generic.cc
View Options

File Metadata

Event Timeline

test_ME_generic.ccNo OneTemporaryActions

test_ME_generic.ccView Options

File Metadata

Event Timeline

test_ME_generic.cc
No OneTemporary
Actions

test_ME_generic.cc
View Options