check_res.cc
No OneTemporary
Actions

Size

5 KB

Referenced Files

None

Subscribers

None

check_res.cc
View Options

	/**
	* \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 <iostream>
	#include <iterator>
	#include <memory>
	#include <string>
	#include <utility>

	#include "HEJ/Config.hh"
	#include "HEJ/CrossSectionAccumulator.hh"
	#include "HEJ/Event.hh"
	#include "HEJ/event_types.hh"
	#include "HEJ/EventReweighter.hh"
	#include "HEJ/EWConstants.hh"
	#include "HEJ/Fraction.hh"
	#include "HEJ/HiggsCouplingSettings.hh"
	#include "HEJ/Mixmax.hh"
	#include "HEJ/Parameters.hh"
	#include "HEJ/ScaleFunction.hh"
	#include "HEJ/stream.hh"

	#include "fastjet/JetDefinition.hh"

	#include "LHEF/LHEF.h"

	namespace HEJ { struct RNG; }

	namespace {
	const fastjet::JetDefinition JET_DEF{fastjet::kt_algorithm, 0.4};
	const fastjet::JetDefinition BORN_JET_DEF{JET_DEF};
	constexpr double BORN_JETPTMIN = 30;
	const HEJ::Fraction<double> SOFT_PT_REGULATOR{0.1};
	constexpr double JETPTMIN = 35;
	constexpr bool LOG_CORR = false;
	constexpr std::size_t NUM_TRIES = 100;
	constexpr HEJ::ParticleProperties WPROP{80.385, 2.085};
	constexpr HEJ::ParticleProperties ZPROP{91.187, 2.495};
	constexpr HEJ::ParticleProperties HPROP{125, 0.004165};
	constexpr double VEV = 246.2196508;
	using EventTreatment = HEJ::EventTreatment;
	using namespace HEJ::event_type;
	HEJ::EventTreatMap TREAT{
	{no_2_jets, EventTreatment::discard},
	{bad_final_state, EventTreatment::discard},
	{non_resummable, EventTreatment::discard},
	{unof, EventTreatment::discard},
	{unob, EventTreatment::discard},
	{qqxexb, EventTreatment::discard},
	{qqxexf, EventTreatment::discard},
	{qqxmid, EventTreatment::discard},
	{FKL, EventTreatment::reweight}
	};

	bool correct_colour(HEJ::Event const & ev){
	if(!HEJ::event_type::is_resummable(ev.type()))
	return true;
	return ev.is_leading_colour();
	}
	} // namespace

	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.soft_pt_regulator = SOFT_PT_REGULATOR;
	HEJ::MatrixElementConfig ME_conf;
	ME_conf.log_correction = LOG_CORR;
	ME_conf.Higgs_coupling = HEJ::HiggsCouplingSettings{};
	ME_conf.ew_parameters.set_vevWZH(VEV, WPROP, ZPROP, HPROP);
	HEJ::EventReweighterConfig conf;
	conf.psp_config = std::move(psp_conf);
	conf.ME_config = std::move(ME_conf);
	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.
	};
	std::shared_ptr<HEJ::RNG> ran{std::make_shared<HEJ::Mixmax>()};
	HEJ::EventReweighter hej{reader.heprup, std::move(scale_gen), conf, ran};

	HEJ::CrossSectionAccumulator xs;
	do{
	auto ev_data = HEJ::Event::EventData{reader.hepeup};
	shuffle_particles(ev_data);
	ev_data.reconstruct_intermediate();
	HEJ::Event ev{
	ev_data.cluster(
	BORN_JET_DEF, BORN_JETPTMIN
	)
	};
	auto resummed_events = hej.reweight(ev, NUM_TRIES);
	for(auto const & res_ev: resummed_events) {
	ASSERT(correct_colour(res_ev));
	ASSERT(std::isfinite(res_ev.central().weight));
	// we fill the xs uncorrelated since we only want to test the uncertainty
	// of the resummation
	xs.fill(res_ev);
	}
	} while(reader.readEvent());
	const double xsec = xs.total().value;
	const double xsec_err = std::sqrt(xs.total().error);
	const double significance =
	std::abs(xsec - xsec_ref) / std::sqrt( xsec_errxsec_err + tolerancetolerance );
	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;
	}
	return EXIT_SUCCESS;
	}

File Metadata

Mime Type: text/x-c
Expires: Tue, Sep 30, 6:12 AM (8 h, 44 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 6561635
Default Alt Text: check_res.cc (5 KB)

check_res.ccNo OneTemporaryActions

check_res.ccView Options

File Metadata

Event Timeline

check_res.cc
No OneTemporary
Actions

check_res.cc
View Options