diff --git a/t/check_res.cc b/t/check_res.cc index 895be3a..a8a2546 100644 --- a/t/check_res.cc +++ b/t/check_res.cc @@ -1,163 +1,165 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019 * \copyright GPLv2 or later */ #include <iostream> +#include <math.h> #include "LHEF/LHEF.h" #include "HEJ/Event.hh" #include "HEJ/EventReweighter.hh" #include "HEJ/Mixmax.hh" #include "HEJ/stream.hh" #define ASSERT(x) if(!(x)) { \ std::cerr << "Assertion '" #x "' failed.\n"; \ return EXIT_FAILURE; \ } namespace{ const fastjet::JetDefinition jet_def{fastjet::kt_algorithm, 0.4}; const fastjet::JetDefinition Born_jet_def{jet_def}; constexpr double Born_jetptmin = 30; constexpr double extpartonptmin = 30; constexpr double max_ext_soft_pt_fraction = 0.1; constexpr double jetptmin = 35; constexpr bool log_corr = false; using EventTreatment = HEJ::EventTreatment; using namespace HEJ::event_type; HEJ::EventTreatMap treat{ {no_2_jets, EventTreatment::discard}, {bad_final_state, EventTreatment::discard}, {FixedOrder, EventTreatment::discard}, {unof, EventTreatment::discard}, {unob, EventTreatment::discard}, {qqxexb, EventTreatment::discard}, {qqxexf, EventTreatment::discard}, {qqxmid, EventTreatment::discard}, {FKL, EventTreatment::reweight} }; /// true if colour is allowed for particle bool correct_colour(HEJ::Particle const & part){ if(HEJ::is_AWZH_boson(part) && !part.colour) return true; if(!part.colour) return false; int const colour = part.colour->first; int const anti_colour = part.colour->second; if(part.type == HEJ::ParticleID::gluon) return colour != anti_colour && colour > 0 && anti_colour > 0; if(HEJ::is_quark(part)) return anti_colour == 0 && colour > 0; return colour == 0 && anti_colour > 0; } bool correct_colour(HEJ::Event const & ev){ if(!HEJ::event_type::is_HEJ(ev.type())) return true; for(auto const & part: ev.incoming()){ if(!correct_colour(part)) return false; } for(auto const & part: ev.outgoing()){ if(!correct_colour(part)) return false; } return true; } }; int main(int argn, char** argv) { if(argn == 5 && std::string(argv[4]) == "unof"){ --argn; treat[unof] = EventTreatment::reweight; treat[unob] = EventTreatment::discard; treat[FKL] = EventTreatment::discard; } if(argn == 5 && std::string(argv[4]) == "unob"){ --argn; treat[unof] = EventTreatment::discard; treat[unob] = EventTreatment::reweight; treat[FKL] = EventTreatment::discard; } else if(argn == 5 && std::string(argv[4]) == "splitf"){ --argn; treat[qqxexb] = EventTreatment::discard; treat[qqxexf] = EventTreatment::reweight; treat[FKL] = EventTreatment::discard; } else if(argn == 5 && std::string(argv[4]) == "splitb"){ --argn; treat[qqxexb] = EventTreatment::reweight; treat[qqxexf] = EventTreatment::discard; treat[FKL] = EventTreatment::discard; } else if(argn == 5 && std::string(argv[4]) == "qqxmid"){ --argn; treat[qqxmid] = EventTreatment::reweight; treat[FKL] = EventTreatment::discard; } if(argn != 4){ std::cerr << "Usage: check_res eventfile xsection tolerance [uno]"; return EXIT_FAILURE; } const double xsec_ref = std::stod(argv[2]); const double tolerance = std::stod(argv[3]); HEJ::istream in{argv[1]}; LHEF::Reader reader{in}; HEJ::PhaseSpacePointConfig psp_conf; psp_conf.jet_param = HEJ::JetParameters{jet_def, jetptmin}; psp_conf.min_extparton_pt = extpartonptmin; psp_conf.max_ext_soft_pt_fraction = max_ext_soft_pt_fraction; HEJ::MatrixElementConfig ME_conf; ME_conf.log_correction = log_corr; ME_conf.Higgs_coupling = HEJ::HiggsCouplingSettings{}; HEJ::EventReweighterConfig conf; conf.psp_config = std::move(psp_conf); conf.ME_config = std::move(ME_conf); conf.jet_param = psp_conf.jet_param; conf.treat = treat; reader.readEvent(); const bool has_Higgs = std::find( begin(reader.hepeup.IDUP), end(reader.hepeup.IDUP), 25 ) != end(reader.hepeup.IDUP); const double mu = has_Higgs?125.:91.188; HEJ::ScaleGenerator scale_gen{ {{std::to_string(mu), HEJ::FixedScale{mu}}}, {}, 1. }; HEJ::Mixmax ran{}; HEJ::EventReweighter hej{reader.heprup, std::move(scale_gen), conf, ran}; double xsec = 0.; double xsec_err = 0.; do{ auto ev_data = HEJ::Event::EventData{reader.hepeup}; ev_data.reconstruct_intermediate(); HEJ::Event ev{ ev_data.cluster( Born_jet_def, Born_jetptmin ) }; auto resummed_events = hej.reweight(ev, 100); for(auto const & ev: resummed_events) { ASSERT(correct_colour(ev)); + ASSERT(isfinite(ev.central().weight)); xsec += ev.central().weight; xsec_err += ev.central().weight*ev.central().weight; } } while(reader.readEvent()); xsec_err = std::sqrt(xsec_err); const double significance = std::abs(xsec - xsec_ref) / std::sqrt( xsec_err*xsec_err + tolerance*tolerance ); std::cout << xsec_ref << " +/- " << tolerance << " ~ " << xsec << " +- " << xsec_err << " => " << significance << " sigma\n"; if(significance > 3.){ std::cerr << "Cross section is off by over 3 sigma!\n"; return EXIT_FAILURE; } }