No OneTemporary
Actions

Size

38 KB

Subscribers

None

View Options

	diff --git a/src/EventReweighter.cc b/src/EventReweighter.cc
	index ca61aa6..f3563db 100644
	--- a/src/EventReweighter.cc
	+++ b/src/EventReweighter.cc
	@@ -1,326 +1,343 @@
	#include "RHEJ/EventReweighter.hh"

	#include <string>
	#include <unordered_map>

	#include "RHEJ/PhaseSpacePoint.hh"
	#include "RHEJ/PDG_codes.hh"
	#include "RHEJ/utility.hh"

	namespace RHEJ{
	using EventType = event_type::EventType;

	namespace {

	static_assert(
	std::numeric_limits<double>::has_quiet_NaN,
	"no quiet NaN for double"
	);
	constexpr double NaN = std::numeric_limits<double>::quiet_NaN();

	UnclusteredEvent to_UnclusteredEvent(PhaseSpacePoint const & psp){
	UnclusteredEvent result;
	result.incoming = psp.incoming();
	std::sort(
	begin(result.incoming), end(result.incoming),
	[](Particle o1, Particle o2){return o1.p.pz()<o2.p.pz();}
	);
	assert(result.incoming.size() == 2);
	result.outgoing = psp.outgoing();
	assert(
	std::is_sorted(
	begin(result.outgoing), end(result.outgoing),
	rapidity_less{}
	)
	);
	assert(result.outgoing.size() >= 2);
	result.decays = psp.decays();
	result.central.mur = NaN;
	result.central.muf = NaN;
	result.central.weight = psp.weight();
	return result;
	}

	} // namespace anonymous

	EventReweighter::EventReweighter(
	LHEF::HEPRUP const & heprup,
	ScaleGenerator scale_gen,
	EventReweighterConfig conf,
	RHEJ::RNG & ran
	):
	EventReweighter{
	RHEJ::Beam{
	heprup.EBMUP.first,
	{{
	static_cast<RHEJ::ParticleID>(heprup.IDBMUP.first),
	static_cast<RHEJ::ParticleID>(heprup.IDBMUP.second)
	}}
	},
	heprup.PDFSUP.first,
	std::move(scale_gen),
	std::move(conf),
	ran
	}
	{
	if(heprup.EBMUP.second != E_beam_){
	throw std::invalid_argument(
	"asymmetric beam: " + std::to_string(E_beam_)
	+ " ---> <--- " + std::to_string(heprup.EBMUP.second)
	);
	};
	if(heprup.PDFSUP.second != pdf_.id()){
	throw std::invalid_argument(
	"conflicting PDF ids: " + std::to_string(pdf_.id())
	+ " vs. " + std::to_string(heprup.PDFSUP.second)
	);
	}
	}

	EventReweighter::EventReweighter(
	Beam beam,
	int pdf_id,
	ScaleGenerator scale_gen,
	EventReweighterConfig conf,
	RHEJ::RNG & ran
	):
	param_{std::move(conf)},
	E_beam_{beam.E},
	pdf_{pdf_id, beam.type.front(), beam.type.back()},
	MEt2_{
	[this](double mu){ return pdf_.Halphas(mu); },
	param_.ME_config
	},
	scale_gen_(std::move(scale_gen)),
	ran_{ran}
	{}

	PDF const & EventReweighter::pdf() const{
	return pdf_;
	}

	std::vector<Event> EventReweighter::reweight(
	Event const & input_ev, int num_events
	){
	auto res_events = gen_res_events(input_ev, num_events);
	if(res_events.empty()) return {};
	for(auto & event: res_events) event = scale_gen_(event);
	return rescale(input_ev, std::move(res_events));
	}

	/**
	* \brief main generation/reweighting function:
	* generate phase space points and divide out Born factors
	*/
	std::vector<Event> EventReweighter::gen_res_events(
	Event const & ev,
	int phase_space_points
	){
	assert(ev.variations().empty());

	switch(param_.treat.at(ev.type())){
	case EventTreatment::discard: return {};
	case EventTreatment::keep:
	if(! jets_pass_resummation_cuts(ev)) return {};
	else return {ev};
	default:;
	}
	const double Born_shat = shat(ev);

	std::vector<Event> resummation_events;
	for(int psp_number = 0; psp_number < phase_space_points; ++psp_number){
	PhaseSpacePoint psp{ev, param_.psp_config, ran_};
	if(psp.weight() == 0.) continue;
	if(psp.incoming()[0].E() > E_beam_ \|\| psp.incoming()[1].E() > E_beam_) continue;

	resummation_events.emplace_back(
	to_UnclusteredEvent(std::move(psp)),
	param_.jet_param.def, param_.jet_param.min_pt
	);
	auto & new_event = resummation_events.back();
	assert(new_event.variations().empty());
	new_event.central().mur = ev.central().mur;
	new_event.central().muf = ev.central().muf;
	const double resum_shat = shat(new_event);
	new_event.central().weight = ev.central().weightBorn_shat*Born_shat/
	(phase_space_pointsresum_shatresum_shat);
	}
	return resummation_events;
	}

	std::vector<Event> EventReweighter::rescale(
	Event const & Born_ev,
	std::vector<Event> events
	) const{
	const double Born_pdf = pdf_factors(Born_ev).central;
	const double Born_ME = tree_matrix_element(Born_ev);

	for(auto & cur_event: events){
	const auto pdf = pdf_factors(cur_event);
	assert(pdf.variations.size() == cur_event.variations().size());
	const auto ME = matrix_elements(cur_event);
	assert(ME.variations.size() == cur_event.variations().size());
	cur_event.central().weight = pdf.centralME.central/(Born_pdf*Born_ME);
	for(size_t i = 0; i < cur_event.variations().size(); ++i){
	cur_event.variations(i).weight *=
	pdf.variations[i]ME.variations[i]/(Born_pdfBorn_ME);
	}
	}
	return events;
	};

	/**
	* \brief Do the Jets pass the resummation Cuts?
	*
	* @param ev Event in Question
	* @returns 0 or 1 depending on if ev passes Jet Cuts
	*/
	bool EventReweighter::jets_pass_resummation_cuts(
	Event const & ev
	) const{
	const auto out_as_PseudoJet = to_PseudoJet(filter_partons(ev.outgoing()));
	fastjet::ClusterSequence cs{out_as_PseudoJet, param_.jet_param.def};
	return cs.inclusive_jets(param_.jet_param.min_pt).size() == ev.jets().size();
	}



	/**
	* \brief pdf_factors Function
	*
	* @param ev Event in Question
	* @returns EventFactor due to PDFs
	*
	* Calculates the Central value and the variation due
	* to the PDF choice made.
	*/
	EventReweighter::EventFactors
	EventReweighter::pdf_factors(Event const & ev) const{
	auto const & a = ev.incoming().front();
	auto const & b = ev.incoming().back();
	const double xa = a.p.e()/E_beam_;
	const double xb = b.p.e()/E_beam_;

	EventFactors result;
	std::unordered_map<double, double> known_pdf;
	result.central =
	pdf_.pdfpt(0,xa,ev.central().muf,a.type)*
	pdf_.pdfpt(1,xb,ev.central().muf,b.type);
	known_pdf.emplace(ev.central().muf, result.central);

	result.variations.reserve(ev.variations().size());
	for(auto const & ev_param: ev.variations()){
	const double muf = ev_param.muf;
	auto cur_pdf = known_pdf.find(muf);
	if(cur_pdf == known_pdf.end()){
	cur_pdf = known_pdf.emplace(
	muf,
	pdf_.pdfpt(0,xa,muf,a.type)*pdf_.pdfpt(1,xb,muf,b.type)
	).first;
	}
	result.variations.emplace_back(cur_pdf->second);
	}
	assert(result.variations.size() == ev.variations().size());
	return result;
	}



	/**
	* \brief matrix_elements Function
	*
	* @param ev Event in question
	* @returns EventFactor due to MatrixElements
	*
	* Calculates the Central value and the variation due
	* to the Matrix Element.
	*/
	EventReweighter::EventFactors
	EventReweighter::matrix_elements(Event const & ev) const{
	assert(param_.treat.count(ev.type()) > 0);
	if(param_.treat.find(ev.type())->second == EventTreatment::keep){
	return fixed_order_scale_ME(ev);
	}

	// precompute overall kinematic factor
	const double ME_kin = MEt2_.tree_kin(ev.incoming(), ev.outgoing(), true);

	EventFactors result;
	std::unordered_map<double, double> known_ME;
	result.central = MEt2_(
	ev.central().mur,
	ev.incoming(), ev.outgoing(),
	true
	);
	known_ME.emplace(ev.central().mur, result.central);

	result.variations.reserve(ev.variations().size());
	for(auto const & param: ev.variations()){
	const double mur = param.mur;
	auto cur_ME = known_ME.find(mur);
	if(cur_ME == known_ME.end()){
	const double ME = MEt2_.tree_param(
	mur, ev.incoming(), ev.outgoing()
	)ME_kinMEt2_.virtual_corrections(
	mur, ev.incoming(), ev.outgoing()
	);
	cur_ME = known_ME.emplace(mur, ME).first;
	}
	result.variations.emplace_back(cur_ME->second);
	}
	assert(result.variations.size() == ev.variations().size());
	return result;
	}

	/**
	* \brief Computes the tree level matrix element
	*
	* @param ev Event in Question
	* @returns HEJ approximation to Tree level Matrix Element
	*
	* This computes the HEJ approximation to the tree level FO
	* Matrix element which is used within the LO weighting process.
	*/
	double EventReweighter::tree_matrix_element(Event const & ev) const{
	assert(ev.variations().empty());
	assert(param_.treat.count(ev.type()) > 0);
	if(param_.treat.find(ev.type())->second == EventTreatment::keep){
	return fixed_order_scale_ME(ev).central;
	}
	return MEt2_.tree(
	ev.central().mur,
	ev.incoming(), ev.outgoing(),
	false
	);
	}

	/**
	* \brief Scale-dependent part of fixed-order matrix element
	*
	* @param ev Event in question
	* @returns EventFactor scale variation due to FO-ME.
	*
	* This is only called to compute the scale variation for events where
	* we don't do resummation (e.g. non-FKL).
	* Since at tree level the scale dependence is just due to alpha_s,
	* it is enough to return the alpha_s(mur) factors in the matrix element.
	* The rest drops out in the ratio of (output event ME)/(input event ME),
	* so we never have to compute it.
	*/
	EventReweighter::EventFactors
	EventReweighter::fixed_order_scale_ME(Event const & ev) const{
	- const int alpha_s_power = std::count_if(
	- begin(ev.outgoing()), end(ev.outgoing()),
	- [](Particle const & p){ return is_parton(p); }
	- );
	+ int alpha_s_power = 0;
	+ for(auto const & part: ev.outgoing()){
	+ if(is_parton(part))
	+ ++alpha_s_power;
	+ else {
	+ switch(part.type){
	+ case pid::Higgs: {
	+ alpha_s_power += 2;
	+ break;
	+ }
	+ // TODO
	+ case pid::Wp:
	+ case pid::Wm:
	+ case pid::photon:
	+ case pid::Z:
	+ default:
	+ throw std::logic_error("Emission of boson of unsupported type");
	+ }
	+ }
	+ }
	+
	EventFactors result;
	result.central = pow(pdf_.Halphas(ev.central().mur), alpha_s_power);
	for(auto const & var: ev.variations()){
	result.variations.emplace_back(
	pow(pdf_.Halphas(var.mur), alpha_s_power)
	);
	}
	return result;
	}

	-}
	+} // namespace RHEJ
	diff --git a/src/MatrixElement.cc b/src/MatrixElement.cc
	index 61de326..3cb584e 100644
	--- a/src/MatrixElement.cc
	+++ b/src/MatrixElement.cc
	@@ -1,770 +1,786 @@
	#include "RHEJ/MatrixElement.hh"

	#include <CLHEP/Random/Randomize.h>
	#include <CLHEP/Random/RanluxEngine.h>

	#include "RHEJ/Constants.hh"
	#include "RHEJ/currents.hh"
	#include "RHEJ/PDG_codes.hh"
	#include "RHEJ/uno.hh"
	#include "RHEJ/utility.hh"

	namespace RHEJ{
	//cf. last line of eq. (22) in \ref Andersen:2011hs
	double MatrixElement::omega0(
	double alpha_s, double mur,
	fastjet::PseudoJet const & q_j, double lambda
	) const {
	const double result = - alpha_sN_C/M_PIlog(q_j.perp2()/(lambda*lambda));
	if(! param_.log_correction) return result;
	// use alpha_s(sqrt(q_j*lambda)), evolved to mur
	return (
	1. + alpha_s/(4.M_PI)beta0log(murmur/(q_j.perp()*lambda))
	)*result;
	}

	double MatrixElement::virtual_corrections(
	double mur,
	std::array<Particle, 2> const & in,
	std::vector<Particle> const & out
	) const{
	fastjet::PseudoJet const & pa = in.front().p;
	#ifndef NDEBUG
	fastjet::PseudoJet const & pb = in.back().p;
	double const norm = (in.front().p + in.back().p).E();
	#endif

	assert(std::is_sorted(out.begin(), out.end(), rapidity_less{}));
	assert(out.size() >= 2);
	assert(pa.pz() < pb.pz());

	fastjet::PseudoJet q = pa - out[0].p;
	size_t first_idx = 0;
	size_t last_idx = out.size() - 1;
	// if there is a Higgs or unordered gluon outside the extremal partons
	// then it is not part of the FKL ladder and does not contribute
	// to the virtual corrections
	if(out.front().type == pid::Higgs \|\| has_unob_gluon(in, out)){
	q -= out[1].p;
	++first_idx;
	}
	if(out.back().type == pid::Higgs \|\| has_unof_gluon(in, out)){
	--last_idx;
	}

	double exponent = 0;
	const double alpha_s = alpha_s_(mur);
	for(size_t j = first_idx; j < last_idx; ++j){
	exponent += omega0(alpha_s, mur, q, CLAMBDA)*(
	out[j+1].rapidity() - out[j].rapidity()
	);
	q -= out[j+1].p;
	}
	assert(
	nearby(q, -1*pb, norm)
	\|\| out.back().type == pid::Higgs
	\|\| has_unof_gluon(in, out)
	);
	return exp(exponent);
	}
	} // namespace RHEJ

	namespace {
	//! Lipatov vertex for partons emitted into extremal jets
	double C2Lipatov(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv,
	CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2)
	{
	CLHEP::HepLorentzVector temptrans=-(qav+qbv);
	CLHEP::HepLorentzVector p5=qav-qbv;
	CLHEP::HepLorentzVector CL=temptrans
	+ p1(qav.m2()/p5.dot(p1) + 2.p5.dot(p2)/p1.dot(p2))
	- p2(qbv.m2()/p5.dot(p2) + 2.p5.dot(p1)/p1.dot(p2));

	// cout << "#Fadin qa : "<<qav<<endl;
	// cout << "#Fadin qb : "<<qbv<<endl;
	// cout << "#Fadin p1 : "<<p1<<endl;
	// cout << "#Fadin p2 : "<<p2<<endl;
	// cout << "#Fadin p5 : "<<p5<<endl;
	// cout << "#Fadin Gauge Check : "<< CL.dot(p5)<<endl;
	// cout << "#Fadin C2L : "<< -CL.dot(CL)<<" "<<-CL.dot(CL)/(qav.m2()*qbv.m2())/(4./p5.perp2())<<endl;

	// TODO can this dead test go?
	// if (-CL.dot(CL)<0.)
	// if (fabs(CL.dot(p5))>fabs(CL.dot(CL))) // not sufficient!
	// return 0.;
	// else
	return -CL.dot(CL);
	}

	//! Lipatov vertex with soft subtraction for partons emitted into extremal jets
	double C2Lipatovots(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv,
	CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2)
	{
	double kperp=(qav-qbv).perp();
	if (kperp>RHEJ::CLAMBDA)
	return C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2());
	else {
	double Cls=(C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2()));
	return Cls-4./(kperp*kperp);
	}
	}

	//! Lipatov vertex
	double C2Lipatov(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv,
	CLHEP::HepLorentzVector pim, CLHEP::HepLorentzVector pip,
	CLHEP::HepLorentzVector pom, CLHEP::HepLorentzVector pop) // B
	{
	CLHEP::HepLorentzVector temptrans=-(qav+qbv);
	CLHEP::HepLorentzVector p5=qav-qbv;
	CLHEP::HepLorentzVector CL=temptrans
	+ qav.m2()(1./p5.dot(pip)pip + 1./p5.dot(pop)*pop)/2.
	- qbv.m2()(1./p5.dot(pim)pim + 1./p5.dot(pom)*pom)/2.
	+ ( pip*(p5.dot(pim)/pip.dot(pim) + p5.dot(pom)/pip.dot(pom))
	+ pop*(p5.dot(pim)/pop.dot(pim) + p5.dot(pom)/pop.dot(pom))
	- pim*(p5.dot(pip)/pip.dot(pim) + p5.dot(pop)/pop.dot(pim))
	- pom*(p5.dot(pip)/pip.dot(pom) + p5.dot(pop)/pop.dot(pom)) )/2.;

	return -CL.dot(CL);
	}

	//! Lipatov vertex with soft subtraction
	double C2Lipatovots(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv,
	CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pb,
	CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2)
	{
	double kperp=(qav-qbv).perp();
	if (kperp>RHEJ::CLAMBDA)
	return C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2());
	else {
	double Cls=(C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2()));
	double temp=Cls-4./(kperp*kperp);
	return temp;
	}
	}

	/** Matrix element squared for tree-level current-current scattering
	* @param aptype Particle a PDG ID
	* @param bptype Particle b PDG ID
	* @param pn Particle n Momentum
	* @param pb Particle b Momentum
	* @param p1 Particle 1 Momentum
	* @param pa Particle a Momentum
	* @returns ME Squared for Tree-Level Current-Current Scattering
	*/
	double ME_current(
	int aptype, int bptype,
	CLHEP::HepLorentzVector const & pn,
	CLHEP::HepLorentzVector const & pb,
	CLHEP::HepLorentzVector const & p1,
	CLHEP::HepLorentzVector const & pa
	){
	if (aptype==21&&bptype==21) {
	return jM2gg(pn,pb,p1,pa);
	} else if (aptype==21&&bptype!=21) {
	if (bptype > 0)
	return jM2qg(pn,pb,p1,pa);
	else
	return jM2qbarg(pn,pb,p1,pa);
	}
	else if (bptype==21&&aptype!=21) { // ----- \|\| -----
	if (aptype > 0)
	return jM2qg(p1,pa,pn,pb);
	else
	return jM2qbarg(p1,pa,pn,pb);
	}
	else { // they are both quark
	if (bptype>0) {
	if (aptype>0)
	return jM2qQ(pn,pb,p1,pa);
	else
	return jM2qQbar(pn,pb,p1,pa);
	}
	else {
	if (aptype>0)
	return jM2qQbar(p1,pa,pn,pb);
	else
	return jM2qbarQbar(pn,pb,p1,pa);
	}
	}
	throw std::logic_error("unknown particle types");
	}

	/** \brief Matrix element squared for tree-level current-current scattering with Higgs
	* @param aptype Particle a PDG ID
	* @param bptype Particle b PDG ID
	* @param pn Particle n Momentum
	* @param pb Particle b Momentum
	* @param p1 Particle 1 Momentum
	* @param pa Particle a Momentum
	* @param qH t-channel momentum before Higgs
	* @param qHp1 t-channel momentum after Higgs
	* @returns ME Squared for Tree-Level Current-Current Scattering with Higgs
	*/
	double ME_Higgs_current(
	int aptype, int bptype,
	CLHEP::HepLorentzVector const & pn,
	CLHEP::HepLorentzVector const & pb,
	CLHEP::HepLorentzVector const & p1,
	CLHEP::HepLorentzVector const & pa,
	CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
	CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
	double mt, bool include_bottom, double mb
	){
	if (aptype==21&&bptype==21) // gg initial state
	return MH2gg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
	else if (aptype==21&&bptype!=21) {
	if (bptype > 0)
	return MH2qg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4./9.;
	else
	return MH2qbarg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4./9.;
	}
	else if (bptype==21&&aptype!=21) {
	if (aptype > 0)
	return MH2qg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4./9.;
	else
	return MH2qbarg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4./9.;
	}
	else { // they are both quark
	if (bptype>0) {
	if (aptype>0)
	return MH2qQ(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)4.4./(9.*9.);
	else
	return MH2qQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)4.4./(9.*9.);
	}
	else {
	if (aptype>0)
	return MH2qQbar(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)4.4./(9.*9.);
	else
	return MH2qbarQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)4.4./(9.*9.);
	}
	}
	throw std::logic_error("unknown particle types");
	}

	/** \brief Current matrix element squared with Higgs and unordered forward emission
	* @param aptype Particle A PDG ID
	* @param bptype Particle B PDG ID
	* @param punof Unordered Particle Momentum
	* @param pn Particle n Momentum
	* @param pb Particle b Momentum
	* @param p1 Particle 1 Momentum
	* @param pa Particle a Momentum
	* @param qH t-channel momentum before Higgs
	* @param qHp1 t-channel momentum after Higgs
	* @returns ME Squared with Higgs and unordered forward emission
	*/
	double ME_Higgs_current_unof(
	int aptype, int bptype,
	CLHEP::HepLorentzVector const & punof,
	CLHEP::HepLorentzVector const & pn,
	CLHEP::HepLorentzVector const & pb,
	CLHEP::HepLorentzVector const & p1,
	CLHEP::HepLorentzVector const & pa,
	CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
	CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
	double mt, bool include_bottom, double mb
	){
	if (aptype==21&&bptype!=21) {
	if (bptype > 0)
	return jM2unogqHg(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
	else
	return jM2unogqbarHg(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
	}
	else { // they are both quark
	if (bptype>0) {
	if (aptype>0)
	return jM2unogqHQ(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
	else
	return jM2unogqHQbar(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
	}
	else {
	if (aptype>0)
	return jM2unogqbarHQ(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
	else
	return jM2unogqbarHQbar(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
	}
	}
	throw std::logic_error("unknown particle types");
	}

	/** \brief Current matrix element squared with Higgs and unordered backward emission
	* @param aptype Particle A PDG ID
	* @param bptype Particle B PDG ID
	* @param pn Particle n Momentum
	* @param pb Particle b Momentum
	* @param punob Unordered back Particle Momentum
	* @param p1 Particle 1 Momentum
	* @param pa Particle a Momentum
	* @param qH t-channel momentum before Higgs
	* @param qHp1 t-channel momentum after Higgs
	* @returns ME Squared with Higgs and unordered backward emission
	*/
	double ME_Higgs_current_unob(
	int aptype, int bptype,
	CLHEP::HepLorentzVector const & pn,
	CLHEP::HepLorentzVector const & pb,
	CLHEP::HepLorentzVector const & punob,
	CLHEP::HepLorentzVector const & p1,
	CLHEP::HepLorentzVector const & pa,
	CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
	CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
	double mt, bool include_bottom, double mb
	){
	if (bptype==21&&aptype!=21) {
	if (aptype > 0)
	return jM2unobgHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
	else
	return jM2unobgHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
	}
	else { // they are both quark
	if (aptype>0) {
	if (bptype>0)
	return jM2unobqHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
	else
	return jM2unobqbarHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
	}
	else {
	if (bptype>0)
	return jM2unobqHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
	else
	return jM2unobqbarHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
	}
	}
	throw std::logic_error("unknown particle types");
	}

	CLHEP::HepLorentzVector to_HepLorentzVector(RHEJ::Particle const & particle){
	return {particle.p.px(), particle.p.py(), particle.p.pz(), particle.p.E()};
	}
	} // namespace anonymous

	namespace RHEJ{
	MatrixElement::MatrixElement(
	std::function<double (double)> alpha_s,
	MatrixElementConfig conf
	):
	alpha_s_{std::move(alpha_s)},
	param_{std::move(conf)}
	{}

	double MatrixElement::operator()(
	double mur,
	std::array<Particle, 2> const & incoming,
	std::vector<Particle> const & outgoing,
	bool check_momenta
	) const {
	return tree(
	mur,
	incoming, outgoing,
	check_momenta
	)*virtual_corrections(
	mur,
	incoming, outgoing
	);
	}

	double MatrixElement::tree_kin(
	std::array<Particle, 2> const & incoming,
	std::vector<Particle> const & outgoing,
	bool check_momenta
	) const {
	assert(
	std::is_sorted(
	incoming.begin(), incoming.end(),
	[](Particle o1, Particle o2){return o1.p.pz()<o2.p.pz();}
	)
	);
	assert(std::is_sorted(outgoing.begin(), outgoing.end(), rapidity_less{}));

	auto AWZH_boson = std::find_if(
	begin(outgoing), end(outgoing),
	[](Particle const & p){return is_AWZH_boson(p);}
	);

	if(AWZH_boson == end(outgoing)){
	return tree_kin_jets(incoming, outgoing, check_momenta);
	}

	switch(AWZH_boson->type){
	case pid::Higgs: {
	- static constexpr double mH = 125.;
	- const double alpha_s_mH = alpha_s_(mH);
	- return alpha_s_mHalpha_s_mH/(256.pow(M_PI, 5))*tree_kin_Higgs(
	- incoming, outgoing, check_momenta
	- );
	+ return tree_kin_Higgs(incoming, outgoing, check_momenta);
	}
	// TODO
	case pid::Wp:
	case pid::Wm:
	case pid::photon:
	case pid::Z:
	default:
	throw std::logic_error("Emission of boson of unsupported type");
	}
	}

	namespace{
	constexpr int extremal_jet_idx = 1;
	constexpr int no_extremal_jet_idx = 0;

	bool treat_as_extremal(Particle const & parton){
	return parton.p.user_index() == extremal_jet_idx;
	}

	template<class InputIterator>
	double FKL_ladder_weight(
	InputIterator begin_gluon, InputIterator end_gluon,
	CLHEP::HepLorentzVector const & q0,
	CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
	CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pn
	){
	double wt = 1;
	auto qi = q0;
	for(auto gluon_it = begin_gluon; gluon_it != end_gluon; ++gluon_it){
	assert(gluon_it->type == pid::gluon);
	const auto g = to_HepLorentzVector(*gluon_it);
	const auto qip1 = qi - g;

	if(treat_as_extremal(*gluon_it)){
	wt = C2Lipatovots(qip1, qi, pa, pb)C_A;
	} else{
	wt = C2Lipatovots(qip1, qi, pa, pb, p1, pn)C_A;
	}

	qi = qip1;
	}
	return wt;
	}

	} // namespace anonymous

	std::vector<Particle> MatrixElement::tag_extremal_jet_partons(
	std::array<Particle, 2> const & incoming,
	std::vector<Particle> out_partons, bool check_momenta
	) const{
	if(!check_momenta){
	for(auto & parton: out_partons){
	parton.p.set_user_index(no_extremal_jet_idx);
	}
	return out_partons;
	}
	fastjet::ClusterSequence cs(to_PseudoJet(out_partons), param_.jet_param.def);
	const auto jets = sorted_by_rapidity(cs.inclusive_jets(param_.jet_param.min_pt));
	assert(jets.size() >= 2);
	auto most_backward = begin(jets);
	auto most_forward = end(jets) - 1;
	// skip jets caused by unordered emission
	if(has_unob_gluon(incoming, out_partons)){
	assert(jets.size() >= 3);
	++most_backward;
	}
	else if(has_unof_gluon(incoming, out_partons)){
	assert(jets.size() >= 3);
	--most_forward;
	}
	const auto extremal_jet_indices = cs.particle_jet_indices(
	{most_backward, most_forward}
	);
	assert(extremal_jet_indices.size() == out_partons.size());
	for(size_t i = 0; i < out_partons.size(); ++i){
	assert(RHEJ::is_parton(out_partons[i]));
	const int idx = (extremal_jet_indices[i]>=0)?
	extremal_jet_idx:
	no_extremal_jet_idx;
	out_partons[i].p.set_user_index(idx);
	}
	return out_partons;
	}

	double MatrixElement::tree_kin_jets(
	std::array<Particle, 2> const & incoming,
	std::vector<Particle> partons,
	bool check_momenta
	) const {
	partons = tag_extremal_jet_partons(incoming, partons, check_momenta);
	if(has_unob_gluon(incoming, partons) \|\| has_unof_gluon(incoming, partons)){
	throw std::logic_error("unordered emission not implemented for pure jets");
	}

	const auto pa = to_HepLorentzVector(incoming[0]);
	const auto pb = to_HepLorentzVector(incoming[1]);

	const auto p1 = to_HepLorentzVector(partons.front());
	const auto pn = to_HepLorentzVector(partons.back());

	return ME_current(
	incoming[0].type, incoming[1].type,
	pn, pb, p1, pa
	)/(4(N_CN_C - 1))*FKL_ladder_weight(
	begin(partons) + 1, end(partons) - 1,
	pa - p1, pa, pb, p1, pn
	);
	}

	double MatrixElement::tree_kin_Higgs(
	std::array<Particle, 2> const & incoming,
	std::vector<Particle> const & outgoing,
	bool check_momenta
	) const {
	if(has_uno_gluon(incoming, outgoing)){
	return tree_kin_Higgs_between(incoming, outgoing, check_momenta);
	}
	if(outgoing.front().type == pid::Higgs){
	return tree_kin_Higgs_first(incoming, outgoing, check_momenta);
	}
	if(outgoing.back().type == pid::Higgs){
	return tree_kin_Higgs_last(incoming, outgoing, check_momenta);
	}
	return tree_kin_Higgs_between(incoming, outgoing, check_momenta);
	}

	double MatrixElement::MH2_forwardH(
	CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
	CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in,
	CLHEP::HepLorentzVector pH,
	double t1, double t2
	) const{
	ignore(p2out, p2in);
	const double shat = p1in.invariantMass2(p2in);
	// gluon case
	#ifdef RHEJ_BUILD_WITH_QCDLOOP
	if(!param_.Higgs_coupling.use_impact_factors){
	return C_A/C_F1./(16M_PIM_PI)t1/t2*MH2gq_outsideH(
	p1out, p1in, p2out, p2in, pH,
	param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom,
	param_.Higgs_coupling.mb
	);
	}
	#endif
	return 9./2.shatshat*(
	C2gHgp(p1in,p1out,pH) + C2gHgm(p1in,p1out,pH)
	)/(t1*t2);
	}

	double MatrixElement::tree_kin_Higgs_first(
	std::array<Particle, 2> const & incoming,
	std::vector<Particle> const & outgoing,
	bool check_momenta
	) const {
	assert(outgoing.front().type == pid::Higgs);
	if(outgoing[1].type != pid::gluon) {
	assert(incoming.front().type == outgoing[1].type);
	return tree_kin_Higgs_between(incoming, outgoing, check_momenta);
	}
	const auto pH = to_HepLorentzVector(outgoing.front());
	const auto partons = tag_extremal_jet_partons(
	incoming,
	std::vector<Particle>(begin(outgoing) + 1, end(outgoing)),
	check_momenta
	);

	const auto pa = to_HepLorentzVector(incoming[0]);
	const auto pb = to_HepLorentzVector(incoming[1]);

	const auto p1 = to_HepLorentzVector(partons.front());
	const auto pn = to_HepLorentzVector(partons.back());

	const auto q0 = pa - p1 - pH;

	const double t1 = q0.m2();
	const double t2 = (pn - pb).m2();

	double wt = MH2_forwardH(
	p1, pa, pn, pb, pH,
	t1, t2
	)*FKL_ladder_weight(
	begin(partons) + 1, end(partons) - 1,
	q0, pa, pb, p1, pn
	);

	for(auto const & inc: incoming){
	if(inc.type != pid::gluon) wt *= C_F/C_A;
	}
	return wt;
	}

	double MatrixElement::tree_kin_Higgs_last(
	std::array<Particle, 2> const & incoming,
	std::vector<Particle> const & outgoing,
	bool check_momenta
	) const {
	assert(outgoing.back().type == pid::Higgs);
	if(outgoing[outgoing.size()-2].type != pid::gluon) {
	assert(incoming.back().type == outgoing[outgoing.size()-2].type);
	return tree_kin_Higgs_between(incoming, outgoing, check_momenta);
	}
	const auto pH = to_HepLorentzVector(outgoing.back());
	const auto partons = tag_extremal_jet_partons(
	incoming,
	std::vector<Particle>(begin(outgoing), end(outgoing) - 1),
	check_momenta
	);

	const auto pa = to_HepLorentzVector(incoming[0]);
	const auto pb = to_HepLorentzVector(incoming[1]);

	auto p1 = to_HepLorentzVector(partons.front());
	const auto pn = to_HepLorentzVector(partons.back());

	auto q0 = pa - p1;

	const double t1 = q0.m2();
	const double t2 = (pn + pH - pb).m2();

	double wt = MH2_forwardH(
	pn, pb, p1, pa, pH,
	t2, t1
	)*FKL_ladder_weight(
	begin(partons) + 1, end(partons) - 1,
	q0, pa, pb, p1, pn
	);

	for(auto const & inc: incoming){
	if(inc.type != pid::gluon) wt *= C_F/C_A;
	}
	return wt;
	}


	double MatrixElement::tree_kin_Higgs_between(
	std::array<Particle, 2> const & incoming,
	std::vector<Particle> const & outgoing,
	bool check_momenta
	) const {
	const auto the_Higgs = std::find_if(
	begin(outgoing), end(outgoing),
	[](Particle const & s){ return s.type == pid::Higgs; }
	);
	assert(the_Higgs != end(outgoing));
	const auto pH = to_HepLorentzVector(*the_Higgs);
	std::vector<Particle> partons(begin(outgoing), the_Higgs);
	partons.insert(end(partons), the_Higgs + 1, end(outgoing));
	partons = tag_extremal_jet_partons(incoming, partons, check_momenta);

	const auto pa = to_HepLorentzVector(incoming[0]);
	const auto pb = to_HepLorentzVector(incoming[1]);

	auto p1 = to_HepLorentzVector(
	partons[has_unob_gluon(incoming, outgoing)?1:0]
	);
	auto pn = to_HepLorentzVector(
	partons[partons.size() - (has_unof_gluon(incoming, outgoing)?2:1)]
	);

	auto first_after_Higgs = begin(partons) + (the_Higgs-begin(outgoing));
	assert(
	(first_after_Higgs == end(partons) && (
	has_unob_gluon(incoming, outgoing)
	\|\| partons.back().type != pid::gluon
	))
	\|\| first_after_Higgs->rapidity() >= the_Higgs->rapidity()
	);
	assert(
	(first_after_Higgs == begin(partons) && (
	has_unof_gluon(incoming, outgoing)
	\|\| partons.front().type != pid::gluon
	))
	\|\| (first_after_Higgs-1)->rapidity() <= the_Higgs->rapidity()
	);
	// always treat the Higgs as if it were in between the extremal FKL partons
	if(first_after_Higgs == begin(partons)) ++first_after_Higgs;
	else if(first_after_Higgs == end(partons)) --first_after_Higgs;

	// t-channel momentum before Higgs
	auto qH = pa;
	for(auto parton_it = begin(partons); parton_it != first_after_Higgs; ++parton_it){
	qH -= to_HepLorentzVector(*parton_it);
	}

	auto q0 = pa - p1;
	auto begin_ladder = begin(partons) + 1;
	auto end_ladder = end(partons) - 1;

	double current_factor;
	if(has_unob_gluon(incoming, outgoing)){
	current_factor = 9./2.*ME_Higgs_current_unob(
	incoming[0].type, incoming[1].type,
	pn, pb, to_HepLorentzVector(partons.front()), p1, pa, qH, qH - pH,
	param_.Higgs_coupling.mt,
	param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb
	);
	const auto p_unob = to_HepLorentzVector(partons.front());
	q0 -= p_unob;
	p1 += p_unob;
	++begin_ladder;
	}
	else if(has_unof_gluon(incoming, outgoing)){
	current_factor = 9./2.*ME_Higgs_current_unof(
	incoming[0].type, incoming[1].type,
	to_HepLorentzVector(partons.back()), pn, pb, p1, pa, qH, qH - pH,
	param_.Higgs_coupling.mt,
	param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb
	);
	pn += to_HepLorentzVector(partons.back());
	--end_ladder;
	}
	else{
	current_factor = ME_Higgs_current(
	incoming[0].type, incoming[1].type,
	pn, pb, p1, pa, qH, qH - pH,
	param_.Higgs_coupling.mt,
	param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb
	);
	}

	const double ladder_factor = FKL_ladder_weight(
	begin_ladder, first_after_Higgs,
	q0, pa, pb, p1, pn
	)*FKL_ladder_weight(
	first_after_Higgs, end_ladder,
	qH - pH, pa, pb, p1, pn
	);
	return current_factor9./8.ladder_factor;
	}

	double MatrixElement::tree_param_partons(
	double alpha_s, double mur,
	std::vector<Particle> const & partons
	) const{
	const double gs2 = 4.M_PIalpha_s;
	double wt = std::pow(gs2, partons.size());
	if(param_.log_correction){
	// use alpha_s(q_perp), evolved to mur
	assert(partons.size() >= 2);
	for(size_t i = 1; i < partons.size()-1; ++i){
	wt = 1 + alpha_s/(2M_PI)beta0log(mur/partons[i].p.perp());
	}
	}
	return wt;
	}

	double MatrixElement::tree_param(
	double mur,
	std::array<Particle, 2> const & incoming,
	std::vector<Particle> const & outgoing
	) const{
	const double alpha_s = alpha_s_(mur);
	+ auto AWZH_boson = std::find_if(
	+ begin(outgoing), end(outgoing),
	+ [](auto const & p){return is_AWZH_boson(p);}
	+ );
	+ double AWZH_coupling = 1.;
	+ if(AWZH_boson != end(outgoing)){
	+ switch(AWZH_boson->type){
	+ case pid::Higgs: {
	+ AWZH_coupling = alpha_salpha_s/(256.pow(M_PI, 5));
	+ break;
	+ }
	+ // TODO
	+ case pid::Wp:
	+ case pid::Wm:
	+ case pid::photon:
	+ case pid::Z:
	+ default:
	+ throw std::logic_error("Emission of boson of unsupported type");
	+ }
	+ }
	if(has_unob_gluon(incoming, outgoing)){
	- return 4M_PIalpha_s*tree_param_partons(
	+ return AWZH_coupling4M_PIalpha_stree_param_partons(
	alpha_s, mur, filter_partons({begin(outgoing) + 1, end(outgoing)})
	);
	}
	if(has_unof_gluon(incoming, outgoing)){
	- return 4M_PIalpha_s*tree_param_partons(
	+ return AWZH_coupling4M_PIalpha_stree_param_partons(
	alpha_s, mur, filter_partons({begin(outgoing), end(outgoing) - 1})
	);
	}
	- return tree_param_partons(alpha_s, mur, filter_partons(outgoing));
	+ return AWZH_coupling*tree_param_partons(alpha_s, mur, filter_partons(outgoing));
	}

	double MatrixElement::tree(
	double mur,
	std::array<Particle, 2> const & incoming,
	std::vector<Particle> const & outgoing,
	bool check_momenta
	) const {
	return tree_param(mur, incoming, outgoing)*tree_kin(
	incoming, outgoing, check_momenta
	);
	}
	} // namespace RHEJ

File Metadata

Mime Type: text/x-diff
Expires: Tue, Jan 21, 1:44 AM (1 d, 14 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4243493
Default Alt Text: (38 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions