diff --git a/t/test_ME_generic.cc b/t/test_ME_generic.cc
index 8c433cc..546b49c 100644
--- a/t/test_ME_generic.cc
+++ b/t/test_ME_generic.cc
@@ -1,213 +1,280 @@
/**
* \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-2020
* \copyright GPLv2 or later
*/
#include "hej_test.hh"
+#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <fstream>
#include <iostream>
#include <string>
#include "HEJ/Config.hh"
#include "HEJ/Event.hh"
#include "HEJ/EventReader.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/MatrixElement.hh"
+#include "HEJ/stream.hh"
#include "HEJ/YAMLreader.hh"
namespace {
constexpr double alpha_s = 0.118;
- constexpr double ep = 1e-9;
+ 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;
}
+ //! p_z -> -p_z
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
if(ev.decays.size() == 1){
auto decay = ev.decays.begin();
if(ev.outgoing[decay->first].type == HEJ::ParticleID::Z_photon_mix){
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;
}
+
+ void to_anti(HEJ::Particle & part){
+ if(is_anyquark(part)
+ || is_anylepton(part)
+ || abs(part.type) == HEJ::ParticleID::Wp)
+ part.type = static_cast<HEJ::ParticleID>(-part.type);
+ }
+
+ //! p -> -p
+ HEJ::Event::EventData reverse_event(HEJ::Event::EventData ev){
+ for(auto & part: ev.incoming){
+ to_anti(part);
+ auto & p{ part.p };
+ p.reset(-p.px(),-p.py(),-p.pz(),p.E());
+ }
+ for(auto & part: ev.outgoing){
+ to_anti(part);
+ auto & p{ part.p };
+ p.reset(-p.px(),-p.py(),-p.pz(),p.E());
+ }
+ for(auto & decay: ev.decays){
+ for(auto & part: decay.second){
+ to_anti(part);
+ 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(9);
wgt_file << std::scientific;
} 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.;
std::size_t idx_max_ratio = 0;
HEJ::Event ev_max_ratio(HEJ::Event::EventData{}.cluster(
config.resummation_jets.def,0
)
);
double av_ratio = 0;
std::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::EventData data_reverse{reverse_event(data)};
+ shuffle_particles(data_reverse);
+
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
)
};
+ HEJ::Event event_reverse{
+ data_reverse.cluster(
+ config.resummation_jets.def,
+ config.resummation_jets.min_pt
+ )
+ };
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;
}
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){
std::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;
}
+ double ME_reverse;
+ if(component == MEComponent::tree){
+ ME_reverse = ME.tree(event_reverse).central;
+ }else{
+ double sum = 0.;
+ for(size_t i=0; i<ME.virtual_corrections(event_reverse).size(); ++i) {
+ sum += ME.virtual_corrections(event_reverse).at(i).central;
+ }
+ ME_reverse = sum;
+ }
+ if(!std::isfinite(ME_reverse)){
+ std::cerr << "Found non-finite ME ("<< ME_reverse <<")\n" << event_reverse << std::endl;
+ return EXIT_FAILURE;
+ }
+ if(std::abs(our_ME/ME_reverse-1.)>ep_mirror){
+ std::size_t precision(std::cout.precision());
+ std::cerr.precision(16);
+ std::cerr<< "p-reversed ME gives different result " << i << "\n"
+ <<our_ME << " vs " << ME_reverse << " => difference: "
+ << std::abs(our_ME/ME_reverse-1.) << "\n" << event
+ << "\nreversed:\n" << event_reverse << std::endl;
+ std::cerr.precision(precision);
+ // TODO activate
+ // 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);
double diff;
if(ref_ME == 0.){
diff = (our_ME != 0.);
} else {
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 ){
std::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 ) {
std::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;
}