diff --git a/t/test_ME_generic.cc b/t/test_ME_generic.cc
index 8b86cb5..108ed81 100644
--- a/t/test_ME_generic.cc
+++ b/t/test_ME_generic.cc
@@ -1,200 +1,212 @@
/**
* \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());
}
}
// If this is a Z event, we also have to change quarks and leptons to their antiparticle
bool Z_event = false;
if(ev.decays.size() == 1){
for(auto const & decay: ev.decays){
if(ev.outgoing[decay.first].type == HEJ::ParticleID::Z_photon_mix){
Z_event = true;
}
}
}
if(Z_event){
for(auto & part: ev.incoming){
if(is_anyquark(part)){
part.type = (HEJ::ParticleID) -part.type;
}
}
for(auto & part: ev.outgoing){
if(is_anyquark(part)){
part.type = (HEJ::ParticleID) -part.type;
}
}
for(auto & decay: ev.decays){
for(auto & part: decay.second){
part.type = (HEJ::ParticleID) -part.type;
}
}
}
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).at(0).central
- ;
+ double our_ME;
+ if(component == MEComponent::tree){
+ our_ME = ME.tree(event).central;
+ }else{
+ double sum = 0.;
+ for(size_t i=0; i<ME.virtual_corrections(event).size(); ++i) {
+ sum += ME.virtual_corrections(event).at(i).central;
+ }
+ our_ME = sum;
+ }
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).at(0).central
- ;
+ double ME_mirror;
+ if(component == MEComponent::tree){
+ ME_mirror = ME.tree(event_mirror).central;
+ }else{
+ double sum = 0.;
+ for(size_t i=0; i<ME.virtual_corrections(event_mirror).size(); ++i) {
+ sum += ME.virtual_corrections(event_mirror).at(i).central;
+ }
+ ME_mirror = sum;
+ }
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;
}