Page MenuHomeHEPForge
Files F11221425

No OneTemporary
Actions

View Options

diff --git a/t/test_ME_generic.cc b/t/test_ME_generic.cc
index 290fc1c..d649f31 100644
--- a/t/test_ME_generic.cc
+++ b/t/test_ME_generic.cc
@@ -1,132 +1,179 @@
/**
* \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 "HEJ/Event.hh"
#include "HEJ/EventReader.hh"
#include "HEJ/MatrixElement.hh"
#include "HEJ/stream.hh"
#include "HEJ/YAMLreader.hh"
constexpr double alpha_s = 0.118;
constexpr double ep = 1e-5;
+constexpr double ep_mirror = 1e-3;
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);
}
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(),-1*p.pz(),p.E());
+ }
+ for(auto & part: ev.outgoing){
+ auto & p{ part.p };
+ p.reset(p.px(),p.py(),-1*p.pz(),p.E());
+ }
+ for(auto & decay: ev.decays){
+ for(auto & part: decay.second){
+ auto & p{ part.p };
+ p.reset(p.px(),p.py(),-1*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(!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(!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
<< 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-diff
Expires
Wed, May 14, 10:26 AM (1 d, 20 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
5111125
Default Alt Text
(5 KB)

Event Timeline