Page MenuHomeHEPForge
Files F8725678

No OneTemporary
Actions

View Options

diff --git a/FixedOrderGen/include/PhaseSpacePoint.hh b/FixedOrderGen/include/PhaseSpacePoint.hh
index 74257c9..bb0deee 100644
--- a/FixedOrderGen/include/PhaseSpacePoint.hh
+++ b/FixedOrderGen/include/PhaseSpacePoint.hh
@@ -1,92 +1,117 @@
/** \file PhaseSpacePoint.hh
* \brief Contains the PhaseSpacePoint Class
*/
#pragma once
#include <vector>
#include <CLHEP/Random/Randomize.h>
#include <CLHEP/Random/RanluxEngine.h>
#include "fastjet/JetDefinition.hh"
#include "RHEJ/utility.hh"
#include "RHEJ/Event.hh"
+#include "RHEJ/PDG_codes.hh"
namespace HEJFOG{
using RHEJ::Sparticle;
//! A point in resummation phase space
class PhaseSpacePoint{
public:
//! Default PhaseSpacePoint Constructor
PhaseSpacePoint() = default;
//! PhaseSpacePoint Constructor
/**
* @param ev Clustered Jet Event
* @param jet_def The Jet Definition Used
* @param jetptmin Minimum Jet Transverse Momentum
* @param extpartonptmin Minimum Parton Transverse Momentum
* @param max_ext_soft_pt Maximum Parton Transverse Momentum?
*/
PhaseSpacePoint(
int nj, fastjet::JetDefinition jet_def,double jetptmin, double rapmax
);
//! Get Weight Function
/**
* @returns Weight of Event
*/
double weight() const{
return weight_;
}
//! Get Incoming Function
/**
* @returns Incoming Particles
*/
std::array<Sparticle, 2> const & incoming() const{
return incoming_;
}
//! Get Outgoing Function
/**
* @returns Outgoing Particles
*/
std::vector<Sparticle> const & outgoing() const{
return outgoing_;
}
static void reset_ranlux(std::string const & init_file);
static void reset_ranlux(char const * init_file);
private:
+ std::vector<fastjet::PseudoJet> gen_LO_partons(
+ int count, double ptmin, double ptmax
+ );
+ Sparticle gen_boson(
+ RHEJ::ParticleID bosonid
+ );
void rescale_rapidities(
std::vector<fastjet::PseudoJet> & partons,
double ymin, double ymax
);
bool jets_ok(
std::vector<fastjet::PseudoJet> const & Born_jets,
std::vector<fastjet::PseudoJet> const & partons
) const;
void reconstruct_incoming(std::array<Sparticle, 2> const & Born_incoming);
double phase_space_normalisation(
int num_Born_jets
) const;
bool momentum_conserved(double ep) const;
+ template<class Particle>
+ Particle const & most_backward_FKL(
+ std::vector<Particle> const & partons
+ ) const;
+ template<class Particle>
+ Particle const & most_forward_FKL(
+ std::vector<Particle> const & partons
+ ) const;
+ template<class Particle>
+ Particle & most_backward_FKL(std::vector<Particle> & partons) const;
+ template<class Particle>
+ Particle & most_forward_FKL(std::vector<Particle> & partons) const;
+ bool extremal_FKL_ok(
+ std::vector<fastjet::PseudoJet> const & partons
+ ) const;
+
+ bool unob_, unof_;
double weight_;
double jetptmin_;
+
fastjet::JetDefinition jet_def_;
std::array<Sparticle, 2> incoming_;
std::vector<Sparticle> outgoing_;
static CLHEP::Ranlux64Engine ran_;
};
RHEJ::UnclusteredEvent to_UnclusteredEvent(PhaseSpacePoint const & psp);
}
diff --git a/FixedOrderGen/src/PhaseSpacePoint.cc b/FixedOrderGen/src/PhaseSpacePoint.cc
index 25bbc97..cbd99cb 100644
--- a/FixedOrderGen/src/PhaseSpacePoint.cc
+++ b/FixedOrderGen/src/PhaseSpacePoint.cc
@@ -1,474 +1,599 @@
#include "PhaseSpacePoint.hh"
#include <random>
// #include "RHEJ/resummation_jet_momenta.hh"
// #include "RHEJ/Jacobian.hh"
// #include "RHEJ/RNGWrapper.hh"
// #include "RHEJ/uno.hh"
#include "RHEJ/debug.hh"
#include "RHEJ/kinematics.hh"
#include "RHEJ/utility.hh"
+using namespace RHEJ;
+
namespace HEJFOG{
static_assert(
std::numeric_limits<double>::has_quiet_NaN,
"no quiet NaN for double"
);
constexpr double NaN = std::numeric_limits<double>::quiet_NaN();
RHEJ::UnclusteredEvent to_UnclusteredEvent(PhaseSpacePoint const & psp){
RHEJ::UnclusteredEvent result;
result.incoming = psp.incoming();
std::sort(
begin(result.incoming), end(result.incoming),
[](Sparticle o1, Sparticle 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),
RHEJ::rapidity_less{}
)
);
assert(result.outgoing.size() >= 2);
result.central.mur = NaN;
result.central.muf = NaN;
result.central.weight = psp.weight();
return result;
}
namespace{
//generate Ranlux64Engine with fixed, predefined state
/*
* some (all?) of the Ranlux64Engine constructors leave fields
* uninitialised, invoking undefined behaviour. This can be
* circumvented by restoring the state from a file
*/
CLHEP::Ranlux64Engine gen_Ranlux64Engine(){
static const std::string state =
"9876\n"
"0.91280703978419097666\n"
"0.41606065829518357191\n"
"0.99156342622341142601\n"
"0.030922955274050423213\n"
"0.16206278421638486975\n"
"0.76151768001958330956\n"
"0.43765760066092695979\n"
"0.42904698253748563275\n"
"0.11476317525663759511\n"
"0.026620053590963976831\n"
"0.65953715764414511114\n"
"0.30136722624439826745\n"
"3.5527136788005009294e-15 4\n"
"1 202\n";
const std::string file = std::tmpnam(nullptr);
{
std::ofstream out{file};
out << state;
}
CLHEP::Ranlux64Engine result;
result.restoreStatus(file.c_str());
return result;
}
}
CLHEP::Ranlux64Engine PhaseSpacePoint::ran_{gen_Ranlux64Engine()};
void PhaseSpacePoint::reset_ranlux(std::string const & init_file){
reset_ranlux(init_file.c_str());
}
void PhaseSpacePoint::reset_ranlux(char const * init_file){
ran_.restoreStatus(init_file);
}
namespace {
// constexpr int max_jet_user_idx = PhaseSpacePoint::ng_max;
// bool is_nonjet_parton(fastjet::PseudoJet const & parton){
// assert(parton.user_index() != -1);
// return parton.user_index() > max_jet_user_idx;
// }
// bool is_jet_parton(fastjet::PseudoJet const & parton){
// assert(parton.user_index() != -1);
// return parton.user_index() <= max_jet_user_idx;
// }
// user indices for partons with extremal rapidity
constexpr int y_min_user_idx = -3;
constexpr int y_max_user_idx = -2;
}
// namespace{
// double estimate_ng_mean(std::vector<fastjet::PseudoJet> const & Born_jets){
// const double delta_y =
// Born_jets.back().rapidity() - Born_jets.front().rapidity();
// assert(delta_y > 0);
// // Formula derived from fit to 2 jet data
// // for 3 jets: -0.0131111 + (1.39385 + 0.050085*delta_y)*delta_y
// return -0.0213569 + (1.39765 + 0.0498387*delta_y)*delta_y;
// }
// }
// std::vector<fastjet::PseudoJet> PhaseSpacePoint::cluster_jets(
// std::vector<fastjet::PseudoJet> const & partons
// ) const{
// fastjet::ClusterSequence cs(partons, jet_def_);
// return cs.inclusive_jets(jetptmin_);
// }
// bool PhaseSpacePoint::pass_resummation_cuts(
// std::vector<fastjet::PseudoJet> const & jets
// ) const{
// return cluster_jets(jets).size() == jets.size();
// }
// int PhaseSpacePoint::sample_ng(std::vector<fastjet::PseudoJet> const & Born_jets){
// const double ng_mean = estimate_ng_mean(Born_jets);
// std::poisson_distribution<int> dist(ng_mean);
// RNGWrapper<CLHEP::Ranlux64Engine> rng{ran_};
// const int ng = dist(rng);
// assert(ng >= 0);
// assert(ng < ng_max);
// weight_ *= std::tgamma(ng + 1)*std::exp(ng_mean)*std::pow(ng_mean, -ng);
// return ng;
// }
// namespace {
// void copy_AWZH_boson(
// std::vector<Sparticle> const & from,
// std::vector<Sparticle> & to
// ){
// const auto AWZH_boson = std::find_if(
// begin(from), end(from),
// [](Sparticle const & p){ return is_AWZH_boson(p); }
// );
// if(AWZH_boson == end(from)) return;
// auto insertion_point = std::lower_bound(
// begin(to), end(to), *AWZH_boson, rapidity_less{}
// );
// to.insert(insertion_point, *AWZH_boson);
// assert(std::is_sorted(begin(to), end(to), rapidity_less{}));
// }
// }
PhaseSpacePoint::PhaseSpacePoint(
int nj, fastjet::JetDefinition jet_def,double jetptmin, double rapmax
- )
+ ):
+ unob_(false),
+ unof_(false),
+ jetptmin_{jetptmin},
+ jet_def_{jet_def}
{
weight_ = 1;
weight_ /= std::tgamma(nj + 1);
{
outgoing_.reserve(nj + 1); // one slot for possible A, W, Z, H
// use a few variables to check for compilation
{
fastjet::JetDefinition jet_def2(jet_def);
std::cout << " nj: " <<nj<<" jetptmin: "<<jetptmin<<" rapmax : "<<rapmax<<std::endl;
+
+ // sqrt(s)/2 is the maximum pt
+ std::vector<fastjet::PseudoJet> out_partons = gen_LO_partons(nj , jetptmin_,13000./2.);
+
+ if (out_partons.empty())
+ return;
+ // define flavour and incoming states
- // if(&jet_def==NULL)
+ // Make them all gluons
+ for(auto & p: out_partons){
+ outgoing_.emplace_back(Sparticle{pid::gluon, std::move(p)});
+ }
}
- // for(auto & p: out_partons){
- // outgoing_.emplace_back(Sparticle{pid::gluon, std::move(p)});
- // }
+ // call to pdfs and determination of flavour IDs missing here
+ most_backward_FKL(outgoing_).type = pid::gluon;
+ most_forward_FKL(outgoing_).type = pid::gluon;
+
+ // The Higgs
+ auto pH=gen_boson(pid::higgs);
+
+ // overwrite
// most_backward_FKL(outgoing_).type = ev.incoming().front().type;
// most_forward_FKL(outgoing_).type = ev.incoming().back().type;
// copy_AWZH_boson(ev.outgoing(), outgoing_);
// assert(!outgoing_.empty());
- // reconstruct_incoming(ev.incoming());
+ // reconstruct_incoming(ev.incoming());
}
}
+ std::vector<fastjet::PseudoJet> PhaseSpacePoint::gen_LO_partons(
+ int np , double ptmin, double ptmax
+ ){
+ if (np<2) throw std::invalid_argument{"Not enough partons in gen_LO_partons"};
+ // heuristic parameters for pt sampling
+ const double ptpar = ptmin + np/5.;
+ const double temp1 = atan((ptmax - ptmin)/ptpar);
+
+ std::vector<fastjet::PseudoJet> partons(np);
+ for(size_t i = 0; i < (size_t) np; ++i){
+ const double r1 = ran_.flat();
+ const double pt = ptmin + ptpar*tan(r1*temp1);
+ const double temp2 = cos(r1*temp1);
+ const double phi = 2*M_PI*ran_.flat();
+ weight_ *= 2.0*M_PI*pt*ptpar*temp1/(temp2*temp2);
+ // we don't know the allowed rapidity span yet,
+ // set a random value to be rescaled later on
+ const double y = -5.+10.*ran_.flat();
+ weight_*=10.;
+ partons[i].reset_PtYPhiM(pt, y, phi);
+ // Set user index higher than any jet-parton index
+ // in order to assert that these are not inside jets
+ partons[i].set_user_index(i + 1);
+
+ assert(ptmin-1e-5 <= partons[i].pt() && partons[i].pt() <= ptmax+1e-5);
+ }
+
+ // Need to check that at LO, the number of jets = number of partons;
+ // assert(std::all_of(partons.cbegin(), partons.cend(), is_nonjet_parton));
+ fastjet::ClusterSequence cs(partons, jet_def_);
+ auto cluster_jets=cs.inclusive_jets(jetptmin_);
+ if (cluster_jets.size()!=unsigned(np)){
+ weight_=0.0;
+ return {};
+ }
+
+ std::sort(begin(partons), end(partons), rapidity_less{});
+
+ return partons;
+ }
+
+ Sparticle PhaseSpacePoint::gen_boson(
+ RHEJ::ParticleID bosonid
+ ){
+ std::array<double,2> ptrans{0.,0.};
+ for (auto const & parton:outgoing_) {
+ ptrans[0]-=parton.px();
+ ptrans[1]-=parton.py();
+ }
+
+ // The Higgs:
+ // Generate a y Gaussian distributed around 0
+ double lninvr1,r1,r2,temp,a;
+ r1=ran_.flat();
+ r2=ran_.flat();
+ lninvr1=-log(r1);
+ a=1.6;
+ temp=a*sqrt(2.*lninvr1)*cos(2.*M_PI*r2);
+ double y=temp;
+ weight_=weight_*(exp(temp*temp/2./a/a))*sqrt(2.*M_PI)*a;
+
+ // r1=ran.flat();
+ // double sH=flags.mH*(flags.mH + flags.GammaH*tan((M_PI*r1)/2. + (-1. + r1)*atan(flags.mH/flags.GammaH)));
+ double sH=125.*125.;
+
+ double mHperp=sqrt(ptrans[0]*ptrans[0]+ptrans[1]*ptrans[1]+sH);
+ double pz=mHperp*sinh(y);
+ double E=mHperp*cosh(y);
+
+ return Sparticle{bosonid,fastjet::PseudoJet{ptrans[0],ptrans[1],pz,E}};
+ }
+
void PhaseSpacePoint::rescale_rapidities(
std::vector<fastjet::PseudoJet> & partons,
double ymin, double ymax
){
constexpr double ep = 1e-7;
for(auto & parton: partons){
assert(0 <= parton.rapidity() && parton.rapidity() <= 1);
const double dy = ymax - ymin - 2*ep;
const double y = ymin + ep + dy*parton.rapidity();
parton.reset_momentum_PtYPhiM(parton.pt(), y, parton.phi());
weight_ *= dy;
assert(ymin <= parton.rapidity() && parton.rapidity() <= ymax);
}
}
+ template<class Particle>
+ Particle const & PhaseSpacePoint::most_backward_FKL(
+ std::vector<Particle> const & partons
+ ) const{
+ return partons[0 + unob_?1:0];
+ }
+
+ template<class Particle>
+ Particle const & PhaseSpacePoint::most_forward_FKL(
+ std::vector<Particle> const & partons
+ ) const{
+ const size_t idx = partons.size() - 1 - unof_;
+ assert(idx < partons.size());
+ return partons[idx];
+ }
+
+ template<class Particle>
+ Particle & PhaseSpacePoint::most_backward_FKL(
+ std::vector<Particle> & partons
+ ) const{
+ return partons[0 + unob_];
+ }
+
+ template<class Particle>
+ Particle & PhaseSpacePoint::most_forward_FKL(
+ std::vector<Particle> & partons
+ ) const{
+ const size_t idx = partons.size() - 1 - unof_;
+ assert(idx < partons.size());
+ return partons[idx];
+ }
+
// int PhaseSpacePoint::sample_ng_jets(
// int ng, std::vector<fastjet::PseudoJet> const & Born_jets
// ){
// RNGWrapper<CLHEP::Ranlux64Engine> rng{ran_};
// const double p_J = probability_in_jet(Born_jets);
// std::binomial_distribution<> bin_dist(ng, p_J);
// const int ng_J = bin_dist(rng);
// weight_ *= std::pow(p_J, -ng_J)*std::pow(1 - p_J, ng_J - ng);
// return ng_J;
// }
// std::vector<fastjet::PseudoJet>
// PhaseSpacePoint::reshuffle(
// std::vector<fastjet::PseudoJet> const & Born_jets,
// fastjet::PseudoJet const & q
// ){
// if(q == fastjet::PseudoJet{0, 0, 0, 0}) return Born_jets;
// std::vector<fastjet::PseudoJet> jets = resummation_jet_momenta(Born_jets, q);
// if(jets.empty()){
// weight_ = 0;
// return {};
// }
// // transform delta functions to integration over resummation momenta
// weight_ /= Jacobian(jets, q);
// return jets;
// }
// std::vector<int> PhaseSpacePoint::distribute_jet_partons(
// int ng_jets, std::vector<fastjet::PseudoJet> const & jets
// ){
// size_t first_valid_jet = 0;
// size_t num_valid_jets = jets.size();
// const double R_eff = 5./3.*jet_def_.R();
// // if there is an unordered jet too far away from the FKL jets
// // then extra gluon constituents of the unordered jet would
// // violate the FKL rapidity ordering
// if(unob_ && jets[0].delta_R(jets[1]) > R_eff){
// ++first_valid_jet;
// --num_valid_jets;
// }
// else if(unof_ && jets[jets.size()-1].delta_R(jets[jets.size()-2]) > R_eff){
// --num_valid_jets;
// }
// std::vector<int> np(jets.size(), 1);
// for(int i = 0; i < ng_jets; ++i){
// ++np[first_valid_jet + ran_.flat() * num_valid_jets];
// }
// weight_ *= std::pow(num_valid_jets, ng_jets);
// return np;
// }
// std::vector<fastjet::PseudoJet> PhaseSpacePoint::split(
// std::vector<fastjet::PseudoJet> const & jets,
// int ng_jets
// ){
// return split(jets, distribute_jet_partons(ng_jets, jets));
// }
// bool PhaseSpacePoint::pass_extremal_cuts(
// fastjet::PseudoJet const & ext_parton,
// fastjet::PseudoJet const & jet
// ) const{
// if(ext_parton.pt() < extpartonptmin_) return false;
// return (ext_parton - jet).pt()/ext_parton.pt() < max_ext_soft_pt_fraction_;
// }
// std::vector<fastjet::PseudoJet> PhaseSpacePoint::split(
// std::vector<fastjet::PseudoJet> const & jets,
// std::vector<int> const & np
// ){
// assert(! jets.empty());
// assert(jets.size() == np.size());
// assert(pass_resummation_cuts(jets));
// const size_t most_backward_FKL_idx = 0 + unob_;
// const size_t most_forward_FKL_idx = jets.size() - 1 - unof_;
// std::vector<fastjet::PseudoJet> jet_partons;
// // randomly distribute jet gluons among jets
// for(size_t i = 0; i < jets.size(); ++i){
// weight_ *= splitter_.Split(jets[i], np[i]);
// if(weight_ == 0) return {};
// assert(
// std::all_of(
// begin(splitter_.get_jcons()), end(splitter_.get_jcons()),
// is_jet_parton
// )
// );
// const auto first_new_parton = jet_partons.insert(
// end(jet_partons),
// begin(splitter_.get_jcons()), end(splitter_.get_jcons())
// );
// auto extremal = end(jet_partons);
// if((unob_ && i == 0) || i == most_backward_FKL_idx){
// // unordered or FKL backward emission
// extremal = std::min_element(
// first_new_parton, end(jet_partons), rapidity_less{}
// );
// if(i == most_backward_FKL_idx) extremal->set_user_index(y_min_user_idx);
// }
// else if((unof_ && i == jets.size() - 1) || i == most_forward_FKL_idx){
// // unordered or FKL forward emission
// extremal = std::max_element(
// first_new_parton, end(jet_partons), rapidity_less{}
// );
// if(i == most_forward_FKL_idx) extremal->set_user_index(y_max_user_idx);
// }
// if(
// extremal != end(jet_partons)
// && !pass_extremal_cuts(*extremal, jets[i])
// ){
// weight_ = 0;
// return {};
// }
// }
// assert(tagged_FKL_extremal(jet_partons));
// std::sort(begin(jet_partons), end(jet_partons), rapidity_less{});
// if(
// !extremal_FKL_ok(jet_partons)
// || !split_preserved_jets(jets, jet_partons)
// ){
// weight_ = 0.;
// return {};
// }
// return jet_partons;
// }
// bool PhaseSpacePoint::extremal_FKL_ok(
// std::vector<fastjet::PseudoJet> const & partons
// ) const{
// assert(std::is_sorted(begin(partons), end(partons), rapidity_less{}));
// return
// most_backward_FKL(partons).user_index() == y_min_user_idx
// && most_forward_FKL(partons).user_index() == y_max_user_idx;
// }
// bool PhaseSpacePoint::split_preserved_jets(
// std::vector<fastjet::PseudoJet> const & jets,
// std::vector<fastjet::PseudoJet> const & jet_partons
// ) const{
// assert(std::is_sorted(begin(jets), end(jets), rapidity_less{}));
// const auto split_jets = sorted_by_rapidity(cluster_jets(jet_partons));
// // this can happen if two overlapping jets
// // are both split into more than one parton
// if(split_jets.size() != jets.size()) return false;
// for(size_t i = 0; i < split_jets.size(); ++i){
// // this can happen if there are two overlapping jets
// // and a parton is assigned to the "wrong" jet
// if(!nearby_ep(jets[i].rapidity(), split_jets[i].rapidity(), 1e-2)){
// return false;
// }
// }
// return true;
// }
// template<class Particle>
// Particle const & PhaseSpacePoint::most_backward_FKL(
// std::vector<Particle> const & partons
// ) const{
// return partons[0 + unob_];
// }
// template<class Particle>
// Particle const & PhaseSpacePoint::most_forward_FKL(
// std::vector<Particle> const & partons
// ) const{
// const size_t idx = partons.size() - 1 - unof_;
// assert(idx < partons.size());
// return partons[idx];
// }
// template<class Particle>
// Particle & PhaseSpacePoint::most_backward_FKL(
// std::vector<Particle> & partons
// ) const{
// return partons[0 + unob_];
// }
// template<class Particle>
// Particle & PhaseSpacePoint::most_forward_FKL(
// std::vector<Particle> & partons
// ) const{
// const size_t idx = partons.size() - 1 - unof_;
// assert(idx < partons.size());
// return partons[idx];
// }
// namespace{
// bool contains_idx(
// fastjet::PseudoJet const & jet, fastjet::PseudoJet const & parton
// ){
// auto const & constituents = jet.constituents();
// const int idx = parton.user_index();
// return std::find_if(
// begin(constituents), end(constituents),
// [idx](fastjet::PseudoJet const & con){return con.user_index() == idx;}
// ) != end(constituents);
// }
// }
// /**
// * final jet test:
// * - number of jets must match Born kinematics
// * - no partons designated as nonjet may end up inside jets
// * - all other outgoing partons *must* end up inside jets
// * - the extremal (in rapidity) partons must be inside the extremal jets
// * - rapidities must be the same (by construction)
// */
// bool PhaseSpacePoint::jets_ok(
// std::vector<fastjet::PseudoJet> const & Born_jets,
// std::vector<fastjet::PseudoJet> const & partons
// ) const{
// fastjet::ClusterSequence cs(partons, jet_def_);
// const auto jets = sorted_by_rapidity(cs.inclusive_jets(jetptmin_));
// if(jets.size() != Born_jets.size()) return false;
// int in_jet = 0;
// for(size_t i = 0; i < jets.size(); ++i){
// assert(jets[i].has_constituents());
// for(auto && parton: jets[i].constituents()){
// if(is_nonjet_parton(parton)) return false;
// }
// in_jet += jets[i].constituents().size();
// }
// const int expect_in_jet = std::count_if(
// partons.cbegin(), partons.cend(), is_jet_parton
// );
// if(in_jet != expect_in_jet) return false;
// // note that PseudoJet::contains does not work here
// if(! (
// contains_idx(most_backward_FKL(jets), most_backward_FKL(partons))
// && contains_idx(most_forward_FKL(jets), most_forward_FKL(partons))
// )) return false;
// for(size_t i = 0; i < jets.size(); ++i){
// assert(nearby_ep(jets[i].rapidity(), Born_jets[i].rapidity(), 1e-2));
// }
// return true;
// }
namespace{
}
void PhaseSpacePoint::reconstruct_incoming(
std::array<Sparticle, 2> const & Born_incoming
){
std::tie(incoming_[0].p, incoming_[1].p) = incoming_momenta(outgoing_);
for(size_t i = 0; i < incoming_.size(); ++i){
incoming_[i].type = Born_incoming[i].type;
}
assert(momentum_conserved(1e-7));
}
double PhaseSpacePoint::phase_space_normalisation(
int num_Born_jets
) const{
return pow(16*pow(M_PI,3), num_Born_jets);
}
bool PhaseSpacePoint::momentum_conserved(double ep) const{
fastjet::PseudoJet diff;
for(auto const & in: incoming()) diff += in.p;
for(auto const & out: outgoing()) diff -= out.p;
return nearby_ep(diff, fastjet::PseudoJet{}, ep);
}
}
diff --git a/FixedOrderGen/src/main.cc b/FixedOrderGen/src/main.cc
index 4e382cc..d9f6c2c 100644
--- a/FixedOrderGen/src/main.cc
+++ b/FixedOrderGen/src/main.cc
@@ -1,128 +1,129 @@
/**
* Name: main.cc
* Authors: Jeppe R. Andersen
*/
#include <fstream>
#include <algorithm>
#include <memory>
#include <chrono>
#include <iostream>
#include "yaml-cpp/yaml.h"
#include "LHEF/LHEF.h"
#include "RHEJ/CombinedEventWriter.hh"
#include "RHEJ/get_analysis.hh"
#include "RHEJ/utility.hh"
#include "RHEJ/PDF.hh"
//#include "RHEJ/EventReweighter.hh"
#include "RHEJ/config.hh"
#include "RHEJ/stream.hh"
#include "RHEJ/MatrixElement.hh"
#include "RHEJ/LesHouchesWriter.hh"
#include "PhaseSpacePoint.hh"
RHEJ::Config load_config(char const * filename){
try{
return RHEJ::load_config(filename);
}
catch(std::exception const & exc){
std::cerr << "Error: " << exc.what() << '\n';
std::exit(EXIT_FAILURE);
}
}
std::unique_ptr<RHEJ::Analysis> get_analysis(
YAML::Node const & parameters
){
try{
return RHEJ::get_analysis(parameters);
}
catch(std::exception const & exc){
std::cerr << "Failed to load analysis: " << exc.what() << '\n';
std::exit(EXIT_FAILURE);
}
}
int main(int argn, char** argv) {
std::cout << " __ ___ __ ______ __ __ \n";
std::cout << " / / / (_)___ _/ /_ / ____/___ ___ _________ ___ __ / /__ / /______ \n";
std::cout << " / /_/ / / __ `/ __ \\ / __/ / __ \\/ _ \\/ ___/ __ `/ / / / __ / / _ \\/ __/ ___/ \n";
std::cout << " / __ / / /_/ / / / / / /___/ / / / __/ / / /_/ / /_/ / / /_/ / __/ /_(__ ) \n";
std::cout << " /_/ /_/_/\\__, /_/ /_/ /_____/_/ /_/\\___/_/ \\__, /\\__, / \\____/\\___/\\__/____/ \n";
std::cout << " ____///__/ __ ____ ///__//____/ ______ __ \n";
std::cout << " / ____(_) _____ ____/ / / __ \\_________/ /__ _____ / ____/__ ____ ___ _________ _/ /_____ _____\n";
std::cout << " / /_ / / |/_/ _ \\/ __ / / / / / ___/ __ / _ \\/ ___/ / / __/ _ \\/ __ \\/ _ \\/ ___/ __ `/ __/ __ \\/ ___/\n";
std::cout << " / __/ / /> </ __/ /_/ / / /_/ / / / /_/ / __/ / / /_/ / __/ / / / __/ / / /_/ / /_/ /_/ / / \n";
std::cout << " /_/ /_/_/|_|\\___/\\__,_/ \\____/_/ \\__,_/\\___/_/ \\____/\\___/_/ /_/\\___/_/ \\__,_/\\__/\\____/_/ \n";
using clock = std::chrono::system_clock;
if (argn < 2) {
std::cerr << "\n# Usage:\n.FOgen config_file\n";
return EXIT_FAILURE;
}
const auto start_time = clock::now();
// read configuration
const RHEJ::Config config = load_config(argv[1]);
// RHEJ::istream in{argv[2]};
// std::unique_ptr<RHEJ::Analysis> analysis = get_analysis(
// config.analysis_parameters
// );
// assert(analysis != nullptr);
// Generate a matrix element:
RHEJ::MatrixElement ME(config.fixed_order_jets.def,config.fixed_order_jets.min_pt,config.log_correction,config.Higgs_coupling);
RHEJ::PDF pdf{230000,RHEJ::ParticleID::proton,RHEJ::ParticleID::proton};
RHEJ::LesHouchesWriter LHOfile{"HEJ.lh",LHEF::HEPRUP{}};
// Perform N trial generations
int nevent = 0;
while(nevent< config.trials){
if (int(10000.*nevent/config.trials) % 100 == 0)
std::cout << ".";
std::cout.flush();
// Generate phase space point
HEJFOG::PhaseSpacePoint psp{2,config.fixed_order_jets.def,config.fixed_order_jets.min_pt,5.0};
std::cout<<psp.weight()<<std::endl;
+
RHEJ::Event ev{to_UnclusteredEvent(std::move(psp)), config.fixed_order_jets.def, config.fixed_order_jets.min_pt};
// evaluate matrix element on this point
ME.tree(pdf.Halphas(ev.central().mur), ev.central().mur,
ev.incoming(), ev.outgoing(),false);
LHOfile.write(ev);
// RHEJ::Event FO_event{
// RHEJ::UnclusteredEvent{reader.hepeup},
// config.fixed_order_jets.def, config.fixed_order_jets.min_pt,
// };
// auto resummed_events = rhej.reweight(
// FO_event,
// config.trials, config.scale_gen
// );
// const auto ev_class = rhej.last_event_class();
// ++nevent_class[ev_class];
// if(resummed_events.empty()) ++nfailed_class[ev_class];
// for(auto const & ev: resummed_events){
// analysis->fill(ev);
// writer.write(ev);
// }
nevent++;
} // main event loop
std::cout << std::endl;
std::chrono::duration<double> run_time = (clock::now() - start_time);
std::cout << "\nTask Runtime: " << run_time.count() << " seconds.\n";
return 0;
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Jan 21, 2:17 AM (1 d, 21 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4222738
Default Alt Text
(28 KB)

Event Timeline