diff --git a/FixedOrderGen/include/PhaseSpacePoint.hh b/FixedOrderGen/include/PhaseSpacePoint.hh index c3a7bea..2635bf8 100644 --- a/FixedOrderGen/include/PhaseSpacePoint.hh +++ b/FixedOrderGen/include/PhaseSpacePoint.hh @@ -1,187 +1,187 @@ /** \file PhaseSpacePoint.hh * \brief Contains the PhaseSpacePoint Class */ #pragma once #include #include "HEJ/Event.hh" #include "HEJ/Particle.hh" #include "HEJ/PDF.hh" #include "HEJ/PDG_codes.hh" #include "HEJ/RNG.hh" #include "JetParameters.hh" #include "ParticleProperties.hh" #include "Status.hh" namespace HEJFOG{ class Process; using HEJ::Particle; //! A point in resummation phase space class PhaseSpacePoint{ public: //! Default PhaseSpacePoint Constructor PhaseSpacePoint() = default; //! PhaseSpacePoint Constructor /** * @param proc The process to generate * @param jet_properties Jet defintion & cuts * @param pdf The pdf set (used for sampling) * @param E_beam Energie of the beam * @param subl_chance Chance to turn a potentially unordered * emission into an actual one * @param subl_channels Possible subleading channels. * see HEJFOG::Subleading * @param particle_properties Properties of producted boson * * Initially, only FKL phase space points are generated. If the most * extremal emission in any direction is a quark or anti-quark and the * next emission is a gluon, subl_chance is the chance to turn this into * an unordered emission event by swapping the two flavours. At most one * unordered emission will be generated in this way. */ PhaseSpacePoint( Process const & proc, JetParameters const & jet_properties, HEJ::PDF & pdf, double E_beam, double subl_chance, unsigned int subl_channels, ParticlesPropMap const & particles_properties, HEJ::RNG & ran ); //! Get Weight Function /** * @returns Weight of Event */ double weight() const{ return weight_; } Status status() const{ return status_; } //! Get Incoming Function /** * @returns Incoming Particles */ std::array const & incoming() const{ return incoming_; } //! Get Outgoing Function /** * @returns Outgoing Particles */ std::vector const & outgoing() const{ return outgoing_; } std::unordered_map> const & decays() const{ return decays_; } private: /** * \internal * \brief Generate LO parton momentum * * @param count Number of partons to generate * @param is_pure_jets If true ensures momentum conservation in x and y * @param jet_param Jet properties to fulfil * @param max_pt max allowed pt for a parton (typically E_CMS) * @param ran Random Number Generator * * @returns Momentum of partons * * Ensures that each parton is in its own jet. * Generation is independent of parton flavour. Output is sorted in rapidity. */ std::vector gen_LO_partons( int count, bool is_pure_jets, JetParameters const & jet_param, double max_pt, HEJ::RNG & ran ); Particle gen_boson( HEJ::ParticleID bosonid, double mass, double width, HEJ::RNG & ran ); template fastjet::PseudoJet gen_last_momentum( ParticleMomenta const & other_momenta, double mass_square, double y ) const; bool jets_ok( std::vector const & Born_jets, std::vector const & partons ) const; /** * \internal * \brief Generate incoming partons according to the PDF * * @param uf Scale used in the PDF */ void reconstruct_incoming( std::array const & ids, HEJ::PDF & pdf, double E_beam, double uf, HEJ::RNG & ran ); HEJ::ParticleID generate_incoming_id( size_t beam_idx, double x, double uf, HEJ::PDF & pdf, HEJ::RNG & ran ); bool momentum_conserved(double ep) const; HEJ::Particle const & most_backward_FKL( std::vector const & partons ) const; HEJ::Particle const & most_forward_FKL( std::vector const & partons ) const; HEJ::Particle & most_backward_FKL(std::vector & partons) const; HEJ::Particle & most_forward_FKL(std::vector & partons) const; bool extremal_FKL_ok( std::vector const & partons ) const; double random_normal(double stddev, HEJ::RNG & ran); void maybe_turn_to_uno(double chance, HEJ::RNG & ran); std::vector decay_boson( HEJ::Particle const & parent, std::vector const & decays, HEJ::RNG & ran ); Decay select_decay_channel( std::vector const & decays, HEJ::RNG & ran ); double gen_hard_pt( int np, double ptmin, double ptmax, double y, HEJ::RNG & ran ); double gen_soft_pt(int np, double ptmax, HEJ::RNG & ran); double gen_parton_pt( int count, JetParameters const & jet_param, double ptmax, double y, HEJ::RNG & ran ); double weight_; Status status_; std::array incoming_; std::vector outgoing_; //! Particle decays in the format {outgoing index, decay products} std::unordered_map> decays_; }; - HEJ::UnclusteredEvent to_UnclusteredEvent(PhaseSpacePoint const & psp); + HEJ::EventData to_EventData(PhaseSpacePoint const & psp); } diff --git a/FixedOrderGen/src/EventGenerator.cc b/FixedOrderGen/src/EventGenerator.cc index eb99958..4fdcbdf 100644 --- a/FixedOrderGen/src/EventGenerator.cc +++ b/FixedOrderGen/src/EventGenerator.cc @@ -1,80 +1,80 @@ #include "EventGenerator.hh" #include "Process.hh" #include "Beam.hh" #include "JetParameters.hh" #include "PhaseSpacePoint.hh" #include "HEJ/Event.hh" #include "HEJ/config.hh" namespace HEJFOG{ EventGenerator::EventGenerator( Process process, Beam beam, HEJ::ScaleGenerator scale_gen, JetParameters jets, int pdf_id, double subl_change, unsigned int subl_channels, ParticlesPropMap particles_properties, HEJ::HiggsCouplingSettings Higgs_coupling, HEJ::RNG & ran ): pdf_{pdf_id, beam.particles[0], beam.particles[1]}, ME_{ [this](double mu){ return pdf_.Halphas(mu); }, HEJ::MatrixElementConfig{ false, std::move(Higgs_coupling) } }, scale_gen_{std::move(scale_gen)}, process_{std::move(process)}, jets_{std::move(jets)}, beam_{std::move(beam)}, subl_change_{subl_change}, subl_channels_{subl_channels}, particles_properties_{std::move(particles_properties)}, ran_{ran} { } HEJ::Event EventGenerator::gen_event(){ HEJFOG::PhaseSpacePoint psp{ process_, jets_, pdf_, beam_.energy, subl_change_, subl_channels_, particles_properties_, ran_ }; status_ = psp.status(); if(status_ != good) return {}; HEJ::Event ev = scale_gen_( HEJ::Event{ - to_UnclusteredEvent(std::move(psp)), + to_EventData(std::move(psp)), jets_.def, jets_.min_pt } ); const double shat = HEJ::shat(ev); const double xa = (ev.incoming()[0].E()-ev.incoming()[0].pz())/(2.*beam_.energy); const double xb = (ev.incoming()[1].E()+ev.incoming()[1].pz())/(2.*beam_.energy); // evaluate matrix element on this point const auto ME_weights = ME_.tree(ev); ev.central().weight *= ME_weights.central/(shat*shat); ev.central().weight *= pdf_.pdfpt(0,xa,ev.central().muf, ev.incoming()[0].type); ev.central().weight *= pdf_.pdfpt(0,xb,ev.central().muf, ev.incoming()[1].type); for(size_t i = 0; i < ev.variations().size(); ++i){ auto & var = ev.variations(i); var.weight *= ME_weights.variations[i]/(shat*shat); var.weight *= pdf_.pdfpt(0,xa,var.muf, ev.incoming()[0].type); var.weight *= pdf_.pdfpt(0,xb,var.muf, ev.incoming()[1].type); } return ev; } } diff --git a/FixedOrderGen/src/PhaseSpacePoint.cc b/FixedOrderGen/src/PhaseSpacePoint.cc index 4ee90cd..961f6aa 100644 --- a/FixedOrderGen/src/PhaseSpacePoint.cc +++ b/FixedOrderGen/src/PhaseSpacePoint.cc @@ -1,461 +1,461 @@ #include "PhaseSpacePoint.hh" #include #include #include #include "HEJ/exceptions.hh" #include "HEJ/kinematics.hh" #include "HEJ/Particle.hh" #include "HEJ/utility.hh" #include "Process.hh" #include "Subleading.hh" using namespace HEJ; namespace HEJFOG{ static_assert( std::numeric_limits::has_quiet_NaN, "no quiet NaN for double" ); constexpr double NaN = std::numeric_limits::quiet_NaN(); - HEJ::UnclusteredEvent to_UnclusteredEvent(PhaseSpacePoint const & psp){ - HEJ::UnclusteredEvent result; - result.incoming = psp.incoming(); + HEJ::EventData to_EventData(PhaseSpacePoint const & psp){ + HEJ::EventData result; + result.incoming() = psp.incoming(); std::sort( - begin(result.incoming), end(result.incoming), + begin(result.incoming()), end(result.incoming()), [](Particle o1, Particle o2){return o1.p.pz()= 2); - result.decays = psp.decays(); - result.central.mur = NaN; - result.central.muf = NaN; - result.central.weight = psp.weight(); + 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{ bool can_swap_to_uno( HEJ::Particle const & p1, HEJ::Particle const & p2 ){ return is_parton(p1) && p1.type != HEJ::pid::gluon && p2.type == HEJ::pid::gluon; } } void PhaseSpacePoint::maybe_turn_to_uno( double chance, HEJ::RNG & ran ){ assert(outgoing_.size() >= 2); const size_t nout = outgoing_.size(); const bool can_be_uno_backward = can_swap_to_uno( outgoing_[0], outgoing_[1] ); const bool can_be_uno_forward = can_swap_to_uno( outgoing_[nout-1], outgoing_[nout-2] ); if(!can_be_uno_backward && !can_be_uno_forward) return; if(ran.flat() < chance){ weight_ /= chance; if(can_be_uno_backward && can_be_uno_forward){ if(ran.flat() < 0.5){ std::swap(outgoing_[0].type, outgoing_[1].type); } else{ std::swap(outgoing_[nout-1].type, outgoing_[nout-2].type); } weight_ *= 2.; } else if(can_be_uno_backward){ std::swap(outgoing_[0].type, outgoing_[1].type); } else{ assert(can_be_uno_forward); std::swap(outgoing_[nout-1].type, outgoing_[nout-2].type); } } else weight_ /= 1 - chance; } template fastjet::PseudoJet PhaseSpacePoint::gen_last_momentum( ParticleMomenta const & other_momenta, const double mass_square, const double y ) const { std::array pt{0.,0.}; for (auto const & p: other_momenta) { pt[0]-= p.px(); pt[1]-= p.py(); } const double mperp = sqrt(pt[0]*pt[0]+pt[1]*pt[1]+mass_square); const double pz=mperp*sinh(y); const double E=mperp*cosh(y); return {pt[0], pt[1], pz, E}; } PhaseSpacePoint::PhaseSpacePoint( Process const & proc, JetParameters const & jet_param, HEJ::PDF & pdf, double E_beam, double subl_change, unsigned int subl_channels, ParticlesPropMap const & particles_properties, HEJ::RNG & ran ) { assert(proc.njets >= 2); if(proc.boson && particles_properties.find(*(proc.boson)) == particles_properties.end()) throw HEJ::missing_option("Boson " +std::to_string(*(proc.boson))+" can't be generated: missing properties"); status_ = good; weight_ = 1; const int nout = proc.njets + (proc.boson?1:0); outgoing_.reserve(nout); const bool is_pure_jets = !proc.boson; auto partons = gen_LO_partons( proc.njets, is_pure_jets, jet_param, E_beam, ran ); for(auto&& p_out: partons) { outgoing_.emplace_back(Particle{pid::gluon, std::move(p_out)}); } if(status_ != good) return; // create boson if(proc.boson){ const auto & boson_prop = particles_properties.at(*proc.boson); auto boson(gen_boson(*proc.boson, boson_prop.mass, boson_prop.width, ran)); const auto pos = std::upper_bound( begin(outgoing_),end(outgoing_),boson,rapidity_less{} ); outgoing_.insert(pos, std::move(boson)); if(! boson_prop.decays.empty()){ const size_t boson_idx = std::distance(begin(outgoing_), pos); decays_.emplace( boson_idx, decay_boson(outgoing_[boson_idx], boson_prop.decays, ran) ); } } // normalisation of momentum-conserving delta function weight_ *= pow(2*M_PI, 4); reconstruct_incoming(proc.incoming, pdf, E_beam, jet_param.min_pt, ran); if(status_ != good) return; // set outgoing states most_backward_FKL(outgoing_).type = incoming_[0].type; most_forward_FKL(outgoing_).type = incoming_[1].type; if(proc.njets > 2 && subl_channels & Subleading::uno) maybe_turn_to_uno(subl_change, ran); } double PhaseSpacePoint::gen_hard_pt( int np , double ptmin, double ptmax, double y, HEJ::RNG & ran ) { // heuristic parameters for pt sampling const double ptpar = ptmin + np/5.; const double arg_small_y = atan((ptmax - ptmin)/ptpar); const double y_cut = 3.; const double r1 = ran.flat(); if(y < y_cut){ const double pt = ptmin + ptpar*tan(r1*arg_small_y); const double temp = cos(r1*arg_small_y); weight_ *= pt*ptpar*arg_small_y/(temp*temp); return pt; } const double ptpar2 = ptpar/(1 + 5*(y-y_cut)); const double temp = 1. - std::exp((ptmin-ptmax)/ptpar2); const double pt = ptmin - ptpar2*std::log(1-r1*temp); weight_ *= pt*ptpar2*temp/(1-r1*temp); return pt; } double PhaseSpacePoint::gen_soft_pt(int np, double max_pt, HEJ::RNG & ran) { constexpr double ptpar = 4.; const double r = ran.flat(); const double pt = max_pt + ptpar/np*std::log(r); weight_ *= pt*ptpar/(np*r); return pt; } double PhaseSpacePoint::gen_parton_pt( int count, JetParameters const & jet_param, double max_pt, double y, HEJ::RNG & ran ) { constexpr double p_small_pt = 0.02; if(! jet_param.peak_pt) { return gen_hard_pt(count, jet_param.min_pt, max_pt, y, ran); } const double r = ran.flat(); if(r > p_small_pt) { weight_ /= 1. - p_small_pt; return gen_hard_pt(count, *jet_param.peak_pt, max_pt, y, ran); } weight_ /= p_small_pt; const double pt = gen_soft_pt(count, *jet_param.peak_pt, ran); if(pt < jet_param.min_pt) { weight_=0.0; status_ = not_enough_jets; return jet_param.min_pt; } return pt; } std::vector PhaseSpacePoint::gen_LO_partons( int np, bool is_pure_jets, JetParameters const & jet_param, double max_pt, HEJ::RNG & ran ){ if (np<2) throw std::invalid_argument{"Not enough partons in gen_LO_partons"}; weight_ /= pow(16.*pow(M_PI,3),np); weight_ /= std::tgamma(np+1); //remove rapidity ordering std::vector partons; partons.reserve(np); for(int i = 0; i < np; ++i){ const double y = -jet_param.max_y + 2*jet_param.max_y*ran.flat(); weight_ *= 2*jet_param.max_y; const bool is_last_parton = i+1 == np; if(is_pure_jets && is_last_parton) { constexpr double parton_mass_sq = 0.; partons.emplace_back(gen_last_momentum(partons, parton_mass_sq, y)); break; } const double phi = 2*M_PI*ran.flat(); weight_ *= 2.0*M_PI; const double pt = gen_parton_pt(np, jet_param, max_pt, y, ran); if(weight_ == 0.0) return {}; partons.emplace_back(fastjet::PtYPhiM(pt, y, phi)); assert(jet_param.min_pt <= partons[i].pt()); assert(partons[i].pt() <= max_pt+1e-5); } // Need to check that at LO, the number of jets = number of partons; fastjet::ClusterSequence cs(partons, jet_param.def); auto cluster_jets=cs.inclusive_jets(jet_param.min_pt); if (cluster_jets.size()!=unsigned(np)){ weight_=0.0; status_ = not_enough_jets; return {}; } std::sort(begin(partons), end(partons), rapidity_less{}); return partons; } Particle PhaseSpacePoint::gen_boson( HEJ::ParticleID bosonid, double mass, double width, HEJ::RNG & ran ){ if(bosonid != HEJ::pid::Higgs) throw HEJ::not_implemented("Production of boson "+std::to_string(bosonid) +" not implemented yet."); // Usual phase space measure weight_ /= 16.*pow(M_PI, 3); // Generate a y Gaussian distributed around 0 // TODO: magic number only for Higgs const double y = random_normal(1.6, ran); const double r1 = ran.flat(); const double sH = mass*( mass + width*tan(M_PI/2.*r1 + (r1-1.)*atan(mass/width)) ); auto p = gen_last_momentum(outgoing_, sH, y); return Particle{bosonid, std::move(p)}; } Particle const & PhaseSpacePoint::most_backward_FKL( std::vector const & partons ) const{ if(!HEJ::is_parton(partons[0])) return partons[1]; return partons[0]; } Particle const & PhaseSpacePoint::most_forward_FKL( std::vector const & partons ) const{ const size_t last_idx = partons.size() - 1; if(!HEJ::is_parton(partons[last_idx])) return partons[last_idx-1]; return partons[last_idx]; } Particle & PhaseSpacePoint::most_backward_FKL( std::vector & partons ) const{ if(!HEJ::is_parton(partons[0])) return partons[1]; return partons[0]; } Particle & PhaseSpacePoint::most_forward_FKL( std::vector & partons ) const{ const size_t last_idx = partons.size() - 1; if(!HEJ::is_parton(partons[last_idx])) return partons[last_idx-1]; return partons[last_idx]; } void PhaseSpacePoint::reconstruct_incoming( std::array const & ids, HEJ::PDF & pdf, double E_beam, double uf, HEJ::RNG & ran ){ std::tie(incoming_[0].p, incoming_[1].p) = incoming_momenta(outgoing_); // calculate xa, xb const double sqrts=2*E_beam; const double xa=(incoming_[0].p.e()-incoming_[0].p.pz())/sqrts; const double xb=(incoming_[1].p.e()+incoming_[1].p.pz())/sqrts; // abort if phase space point is outside of collider energy reach if (xa>1. || xb>1.){ weight_=0; status_ = too_much_energy; return; } // pick pdfs /** @TODO: * ufa, ufb don't correspond to our final scale choice. * The HEJ scale generators currently expect a full event as input, * so fixing this is not completely trivial */ for(size_t i = 0; i < 2; ++i){ if(ids[i] == pid::proton || ids[i] == pid::p_bar){ incoming_[i].type = generate_incoming_id(i, i?xb:xa, uf, pdf, ran); } else { incoming_[i].type = ids[i]; } } assert(momentum_conserved(1e-7)); } namespace { /// particles are ordered, odd = anti ParticleID index_to_pid(size_t i){ if(!i) return pid::gluon; return static_cast(i%2?(i+1)/2:-i/2); } } HEJ::ParticleID PhaseSpacePoint::generate_incoming_id( size_t const beam_idx, double const x, double const uf, HEJ::PDF & pdf, HEJ::RNG & ran ){ std::array pdf_wt; pdf_wt[0] = fabs(pdf.pdfpt(beam_idx,x,uf,pid::gluon)); double pdftot = pdf_wt[0]; for(size_t i = 1; i < pdf_wt.size(); ++i){ pdf_wt[i] = 4./9.*fabs(pdf.pdfpt(beam_idx,x,uf,index_to_pid(i))); pdftot += pdf_wt[i]; } const double r1 = pdftot * ran.flat(); double sum = 0; for(size_t i=0; i < pdf_wt.size(); ++i){ if (r1 < (sum+=pdf_wt[i])){ weight_*= pdftot/pdf_wt[i]; return index_to_pid(i); } } std::cerr << "Error in choosing incoming parton: "< const & decays, HEJ::RNG & ran ){ double br_total = 0.; for(auto const & decay: decays) br_total += decay.branching_ratio; // adjust weight // this is given by (channel branching ratio)/(chance to pick channel) // where (chance to pick channel) = // (channel branching ratio)/(total branching ratio) weight_ *= br_total; const double r1 = br_total*ran.flat(); double br_sum = 0.; for(auto const & decay: decays){ br_sum += decay.branching_ratio; if(r1 < br_sum) return decay; } throw std::logic_error{"unreachable"}; } std::vector PhaseSpacePoint::decay_boson( HEJ::Particle const & parent, std::vector const & decays, HEJ::RNG & ran ){ const auto channel = select_decay_channel(decays, ran); if(channel.products.size() != 2){ throw HEJ::not_implemented{ "only decays into two particles are implemented" }; } std::vector decay_products(channel.products.size()); for(size_t i = 0; i < channel.products.size(); ++i){ decay_products[i].type = channel.products[i]; } // choose polar and azimuth angle in parent rest frame const double E = parent.m()/2; const double theta = 2.*M_PI*ran.flat(); const double cos_phi = 2.*ran.flat()-1.; const double sin_phi = sqrt(1. - cos_phi*cos_phi); // Know 0 < phi < pi const double px = E*cos(theta)*sin_phi; const double py = E*sin(theta)*sin_phi; const double pz = E*cos_phi; decay_products[0].p.reset(px, py, pz, E); decay_products[1].p.reset(-px, -py, -pz, E); for(auto & particle: decay_products) particle.p.boost(parent.p); return decay_products; } } diff --git a/src/EventReweighter.cc b/src/EventReweighter.cc index 88ac45d..c888370 100644 --- a/src/EventReweighter.cc +++ b/src/EventReweighter.cc @@ -1,275 +1,275 @@ /** * \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie * \date 2019 * \copyright GPLv2 or later */ #include "HEJ/EventReweighter.hh" #include #include #include #include #include #include #include #include "fastjet/ClusterSequence.hh" #include "LHEF/LHEF.h" #include "HEJ/Event.hh" #include "HEJ/exceptions.hh" #include "HEJ/Particle.hh" #include "HEJ/PDG_codes.hh" #include "HEJ/PhaseSpacePoint.hh" #include "HEJ/Weights.hh" namespace HEJ{ using EventType = event_type::EventType; namespace { static_assert( std::numeric_limits::has_quiet_NaN, "no quiet NaN for double" ); constexpr double NaN = std::numeric_limits::quiet_NaN(); - UnclusteredEvent to_UnclusteredEvent(PhaseSpacePoint const & psp){ - UnclusteredEvent result; - result.incoming = psp.incoming(); + EventData to_EventData(PhaseSpacePoint const & psp){ + EventData result; + result.incoming() = psp.incoming(); std::sort( - begin(result.incoming), end(result.incoming), + begin(result.incoming()), end(result.incoming()), [](Particle o1, Particle o2){return o1.p.pz()= 2); - result.decays = psp.decays(); - result.central.mur = NaN; - result.central.muf = NaN; - result.central.weight = psp.weight(); + 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, HEJ::RNG & ran ): EventReweighter{ HEJ::Beam{ heprup.EBMUP.first, {{ static_cast(heprup.IDBMUP.first), static_cast(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, HEJ::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 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)); } std::vector 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 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)), + to_EventData(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().weight*Born_shat*Born_shat/ (phase_space_points*resum_shat*resum_shat); } return resummation_events; } std::vector EventReweighter::rescale( Event const & Born_ev, std::vector 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.central*ME.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_pdf*Born_ME); } } return events; }; 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(); } Weights 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_; Weights result; std::unordered_map 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; } Weights 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); } return MEt2_(ev); } 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; } Weights EventReweighter::fixed_order_scale_ME(Event const & ev) const{ 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 not_implemented("Emission of boson of unsupported type"); } } } Weights 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 HEJ diff --git a/src/bin/HEJ.cc b/src/bin/HEJ.cc index 6793c30..dc60d30 100644 --- a/src/bin/HEJ.cc +++ b/src/bin/HEJ.cc @@ -1,188 +1,188 @@ /** * \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie * \date 2019 * \copyright GPLv2 or later */ #include #include #include #include #include #include #include "yaml-cpp/yaml.h" #include "LHEF/LHEF.h" #include "fastjet/ClusterSequence.hh" #include "HEJ/CombinedEventWriter.hh" #include "HEJ/config.hh" #include "HEJ/CrossSectionAccumulator.hh" #include "HEJ/Event.hh" #include "HEJ/EventReweighter.hh" #include "HEJ/get_analysis.hh" #include "HEJ/make_RNG.hh" #include "HEJ/ProgressBar.hh" #include "HEJ/stream.hh" #include "HEJ/Version.hh" #include "HEJ/YAMLreader.hh" int event_number(std::string const & record){ size_t start = record.rfind("Number of Events"); start = record.find_first_of("123456789", start); if(start == std::string::npos) { throw std::invalid_argument("no event number record found"); } const size_t end = record.find_first_not_of("0123456789", start); return std::stoi(record.substr(start, end - start)); } HEJ::Config load_config(char const * filename){ try{ return HEJ::load_config(filename); } catch(std::exception const & exc){ std::cerr << "Error: " << exc.what() << '\n'; std::exit(EXIT_FAILURE); } } std::unique_ptr get_analysis( YAML::Node const & parameters ){ try{ return HEJ::get_analysis(parameters); } catch(std::exception const & exc){ std::cerr << "Failed to load analysis: " << exc.what() << '\n'; std::exit(EXIT_FAILURE); } } // unique_ptr is a workaround: // HEJ::optional is a better fit, but gives spurious errors with g++ 7.3.0 std::unique_ptr> make_progress_bar( std::vector const & xs ) { if(xs.empty()) return {}; const double Born_xs = std::accumulate(begin(xs), end(xs), 0.); return std::make_unique>(std::cout, Born_xs); } std::string time_to_string(const time_t time){ char s[30]; struct tm * p = localtime(&time); strftime(s, 30, "%a %b %d %Y %H:%M:%S", p); return s; } int main(int argn, char** argv) { using clock = std::chrono::system_clock; if (argn < 3) { std::cerr << "\n# Usage:\n."<< argv[0] <<" config_file input_file\n\n"; return EXIT_FAILURE; } const auto start_time = clock::now(); { std::cout << "Starting " << HEJ::Version::package_name_full() << ", revision " << HEJ::Version::revision() << " (" << time_to_string(clock::to_time_t(start_time)) << ")" << std::endl; } fastjet::ClusterSequence::print_banner(); // read configuration const HEJ::Config config = load_config(argv[1]); HEJ::istream in{argv[2]}; std::unique_ptr analysis = get_analysis( config.analysis_parameters ); assert(analysis != nullptr); LHEF::Reader reader{in}; auto heprup = reader.heprup; heprup.generators.emplace_back(LHEF::XMLTag{}); heprup.generators.back().name = HEJ::Version::package_name(); heprup.generators.back().version = HEJ::Version::String(); HEJ::CombinedEventWriter writer{config.output, std::move(heprup)}; double global_reweight = 1.; int max_events = std::numeric_limits::max(); if(argn > 3){ max_events = std::stoi(argv[3]); const int input_events = event_number(reader.headerBlock); global_reweight = input_events/static_cast(max_events); std::cout << "Processing " << max_events << " out of " << input_events << " events\n"; } HEJ::ScaleGenerator scale_gen{ config.scales.base, config.scales.factors, config.scales.max_ratio }; auto ran = HEJ::make_RNG(config.rng.name, config.rng.seed); assert(ran != nullptr); HEJ::EventReweighter hej{ reader.heprup, std::move(scale_gen), to_EventReweighterConfig(config), *ran }; int nevent = 0; std::array nevent_type{0}, nfailed_type{0}; auto progress = make_progress_bar(reader.heprup.XSECUP); HEJ::CrossSectionAccumulator xs; // Loop over the events in the inputfile while(reader.readEvent()){ // reweight events so that the total cross section is conserved reader.hepeup.setWeight(0, global_reweight * reader.hepeup.weight()); if(nevent == max_events) break; ++nevent; // calculate HEJ weight HEJ::Event FO_event{ - HEJ::UnclusteredEvent{reader.hepeup}, + HEJ::EventData{reader.hepeup}, config.fixed_order_jets.def, config.fixed_order_jets.min_pt, }; auto resummed_events = hej.reweight(FO_event, config.trials); ++nevent_type[FO_event.type()]; if(resummed_events.empty()) ++nfailed_type[FO_event.type()]; for(auto const & ev: resummed_events){ //TODO: move pass_cuts to after phase space point generation if(analysis->pass_cuts(ev, FO_event)){ analysis->fill(ev, FO_event); writer.write(ev); xs.fill(ev); } } if(progress) progress->increment(FO_event.central().weight); } // main event loop std::cout << '\n'; analysis->finalise(); using namespace HEJ::event_type; std::cout<< "Events processed: " << nevent << '\n'; for(size_t ev_type = first_type; ev_type <= last_type; ++ev_type){ std::cout << '\t' << names[ev_type] << ": " << nevent_type[ev_type] << ", failed to reconstruct " << nfailed_type[ev_type] << '\n'; } std::cout << '\n' << xs << '\n'; std::chrono::duration run_time = (clock::now() - start_time); std::cout << "Finished " << HEJ::Version::package_name() << " at " << time_to_string(clock::to_time_t(clock::now())) << "\n=> Runtime: " << run_time.count() << " sec (" << nevent/run_time.count() << " Events/sec).\n"; } diff --git a/t/check_res.cc b/t/check_res.cc index 76cab39..108b3ee 100644 --- a/t/check_res.cc +++ b/t/check_res.cc @@ -1,98 +1,98 @@ #include #include "LHEF/LHEF.h" #include "HEJ/Event.hh" #include "HEJ/EventReweighter.hh" #include "HEJ/Mixmax.hh" #include "HEJ/stream.hh" 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 = std::numeric_limits::infinity(); 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}, {nonHEJ, EventTreatment::discard}, {unof, EventTreatment::discard}, {unob, EventTreatment::discard}, {FKL, EventTreatment::reweight} }; }; int main(int argn, char** argv) { if(argn == 5 && std::string(argv[4]) == "uno"){ --argn; treat[unof] = EventTreatment::reweight; treat[unob] = 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{ HEJ::Event ev{ - HEJ::UnclusteredEvent{reader.hepeup}, + HEJ::EventData{reader.hepeup}, Born_jet_def, Born_jetptmin }; auto resummed_events = hej.reweight(ev, 20); for(auto const & ev: resummed_events) { 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; } } diff --git a/t/test_ME_generic.cc b/t/test_ME_generic.cc index 00e19f0..be54a75 100644 --- a/t/test_ME_generic.cc +++ b/t/test_ME_generic.cc @@ -1,104 +1,104 @@ // Generic tester for the ME for a given set of PSP // reference weights and PSP (as LHE file) have to be given as _individual_ files #include #include "LHEF/LHEF.h" #include "HEJ/MatrixElement.hh" #include "HEJ/Event.hh" #include "HEJ/YAMLreader.hh" #include "HEJ/stream.hh" constexpr double alpha_s = 0.118; constexpr double ep = 1e-6; void dump(HEJ::Event const & ev){ { LHEF::Writer writer{std::cout}; std::cout << std::setprecision(6); writer.hepeup = to_HEPEUP(std::move(ev), nullptr); writer.writeEvent(); } std::cout << "Rapidity ordering:\n"; for(const auto & part: ev.outgoing()){ std::cout << std::setw(2) << part.type << ": "<< std::setw(7) << part.rapidity() << std::endl; } } 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); } HEJ::istream in{argv[3]}; LHEF::Reader reader{in}; 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; double av_ratio = 0; size_t i = 0; while(reader.readEvent()){ ++i; HEJ::Event event{ - HEJ::UnclusteredEvent{reader.hepeup}, + HEJ::EventData{reader.hepeup}, config.resummation_jets.def, config.resummation_jets.min_pt }; const double our_ME = ME.tree(event).central; 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::cout.precision(16); std::cout<< "Large difference in PSP " << i << "\nis: "< difference: " << diff << std::endl; std::cout.precision(precision); dump(event); return EXIT_FAILURE; } } } wgt_file.close(); 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; } diff --git a/t/test_classify.cc b/t/test_classify.cc index c7c93d8..e0b58b4 100644 --- a/t/test_classify.cc +++ b/t/test_classify.cc @@ -1,50 +1,50 @@ #include "LHEF/LHEF.h" #include "HEJ/stream.hh" #include "HEJ/event_types.hh" #include "HEJ/Event.hh" namespace{ constexpr double min_jet_pt = 30.; const fastjet::JetDefinition jet_def{fastjet::kt_algorithm, 0.4}; using namespace HEJ::event_type; static const std::vector results{ unob,FKL,FKL,FKL,FKL,FKL,FKL,unob,FKL,unob,FKL,FKL,FKL,unof,FKL,unob,FKL, FKL,unob,unob,FKL,FKL,unob,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,unof, FKL,FKL,unof,FKL,FKL,FKL,FKL,FKL,unof,FKL,FKL,FKL,unof,FKL,FKL,unob,unof, FKL,unof,FKL,unob,FKL,FKL,unob,FKL,unob,unof,unob,unof,FKL,FKL,FKL,FKL,FKL, FKL,FKL,FKL,FKL,FKL,FKL,FKL,unob,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,unob,FKL, FKL,FKL,FKL,unof,FKL,unob,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,unob,FKL, FKL,FKL,FKL,FKL,unob,FKL,unob,unob,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,unof,unob,FKL }; } int main(int argn, char** argv) { if(argn != 2){ std::cerr << "Usage: test_classify eventfile"; return EXIT_FAILURE; } HEJ::istream in{argv[1]}; LHEF::Reader reader{in}; LHEF::Writer writer{std::cerr}; writer.heprup = reader.heprup; for(auto const & expected: results){ reader.readEvent(); const HEJ::Event ev{ - HEJ::UnclusteredEvent{reader.hepeup}, + HEJ::EventData{reader.hepeup}, jet_def, min_jet_pt }; if(ev.type() != expected){ using HEJ::event_type::names; writer.hepeup = reader.hepeup; std::cerr << "wrong classification of event:\n"; writer.writeEvent(); std::cerr << "classified as " << names[ev.type()] << ", is " << names[expected] << '\n'; return EXIT_FAILURE; } } } diff --git a/t/test_descriptions.cc b/t/test_descriptions.cc index d212ad5..f796b02 100644 --- a/t/test_descriptions.cc +++ b/t/test_descriptions.cc @@ -1,62 +1,62 @@ #include #include #include "HEJ/Event.hh" #include "HEJ/EventReweighter.hh" #include "HEJ/ScaleFunction.hh" #define ASSERT(x) if(!(x)) { \ std::cerr << "Assertion '" #x "' failed.\n"; \ return EXIT_FAILURE; \ } int main() { constexpr double mu = 125.; HEJ::ScaleFunction fun{"125", HEJ::FixedScale{mu}}; ASSERT(fun.name() == "125"); HEJ::ScaleGenerator scale_gen{ {std::move(fun)}, {0.5, 1, 2.}, 2.1 }; - HEJ::UnclusteredEvent tmp; - tmp.outgoing.push_back( + HEJ::EventData tmp; + tmp.outgoing().push_back( {HEJ::ParticleID::gluon, fastjet::PtYPhiM(50., -1., 0.3, 0.)} ); - tmp.outgoing.push_back( + tmp.outgoing().push_back( {HEJ::ParticleID::gluon, fastjet::PtYPhiM(30., 1., -0.3, 0.)} ); HEJ::Event ev{ std::move(tmp), fastjet::JetDefinition{fastjet::kt_algorithm, 0.4}, 20. }; auto rescaled = scale_gen(std::move(ev)); ASSERT(rescaled.central().description->scale_name == "125"); for(auto const & var: rescaled.variations()) { ASSERT(var.description->scale_name == "125"); } ASSERT(rescaled.central().description->mur_factor == 1.); ASSERT(rescaled.central().description->muf_factor == 1.); ASSERT(rescaled.variations(0).description->mur_factor == 1.); ASSERT(rescaled.variations(0).description->muf_factor == 1.); ASSERT(rescaled.variations(1).description->mur_factor == 0.5); ASSERT(rescaled.variations(1).description->muf_factor == 0.5); ASSERT(rescaled.variations(2).description->mur_factor == 0.5); ASSERT(rescaled.variations(2).description->muf_factor == 1.); ASSERT(rescaled.variations(3).description->mur_factor == 1.); ASSERT(rescaled.variations(3).description->muf_factor == 0.5); ASSERT(rescaled.variations(4).description->mur_factor == 1.); ASSERT(rescaled.variations(4).description->muf_factor == 2.); ASSERT(rescaled.variations(5).description->mur_factor == 2.); ASSERT(rescaled.variations(5).description->muf_factor == 1.); ASSERT(rescaled.variations(6).description->mur_factor == 2.); ASSERT(rescaled.variations(6).description->muf_factor == 2.); } diff --git a/t/test_psp.cc b/t/test_psp.cc index 0d6f486..35ce3e7 100644 --- a/t/test_psp.cc +++ b/t/test_psp.cc @@ -1,69 +1,69 @@ #include "LHEF/LHEF.h" #include "HEJ/stream.hh" #include "HEJ/config.hh" #include "HEJ/event_types.hh" #include "HEJ/Event.hh" #include "HEJ/PhaseSpacePoint.hh" #include "HEJ/Ranlux64.hh" namespace{ constexpr int max_trials = 100; constexpr double extpartonptmin = 45.; constexpr double max_ext_soft_pt_fraction = std::numeric_limits::infinity(); const fastjet::JetDefinition jet_def{fastjet::kt_algorithm, 0.4}; constexpr double min_jet_pt = 50; }; int main(int argn, char** argv) { if(argn != 2){ std::cerr << "Usage: test_psp eventfile"; return EXIT_FAILURE; } HEJ::istream in{argv[1]}; LHEF::Reader reader{in}; LHEF::Writer writer{std::cerr}; writer.heprup = reader.heprup; HEJ::PhaseSpacePointConfig conf; conf.jet_param = HEJ::JetParameters{jet_def, min_jet_pt}; conf.min_extparton_pt = extpartonptmin; conf.max_ext_soft_pt_fraction = max_ext_soft_pt_fraction; HEJ::Ranlux64 ran{}; while(reader.readEvent()){ const HEJ::Event ev{ - HEJ::UnclusteredEvent{reader.hepeup}, + HEJ::EventData{reader.hepeup}, jet_def, min_jet_pt }; for(int trial = 0; trial < max_trials; ++trial){ HEJ::PhaseSpacePoint psp{ev, conf, ran}; if(psp.weight() != 0){ - HEJ::UnclusteredEvent tmp_ev; - tmp_ev.incoming = psp.incoming(); - tmp_ev.outgoing = psp.outgoing(); - tmp_ev.central = {0,0,0}; + HEJ::EventData tmp_ev; + tmp_ev.incoming() = psp.incoming(); + tmp_ev.outgoing() = psp.outgoing(); + tmp_ev.central() = {0,0,0}; HEJ::Event out_ev{ std::move(tmp_ev), jet_def, min_jet_pt }; if(out_ev.type() != ev.type()){ using HEJ::event_type::names; std::cerr << "Wrong class of phase space point:\n" "original event of class " << names[ev.type()] << ":\n"; writer.hepeup = reader.hepeup; writer.writeEvent(); std::cerr << "Phase space point of class " << names[out_ev.type()] << ":\n"; writer.hepeup = to_HEPEUP(out_ev, &writer.heprup); writer.writeEvent(); return EXIT_FAILURE; } } } } } diff --git a/t/test_scale_arithmetics.cc b/t/test_scale_arithmetics.cc index 015dd6a..484357a 100644 --- a/t/test_scale_arithmetics.cc +++ b/t/test_scale_arithmetics.cc @@ -1,87 +1,87 @@ // Generic tester for the ME for a given set of PSP // reference weights and PSP (as LHE file) have to be given as _individual_ files #include #include "LHEF/LHEF.h" #include "HEJ/EventReweighter.hh" #include "HEJ/make_RNG.hh" #include "HEJ/Event.hh" #include "HEJ/YAMLreader.hh" #include "HEJ/stream.hh" constexpr double alpha_s = 0.118; constexpr double ep = 1e-13; void dump(HEJ::Event const & ev){ { LHEF::Writer writer{std::cout}; std::cout << std::setprecision(6); writer.hepeup = to_HEPEUP(std::move(ev), nullptr); writer.writeEvent(); } std::cout << "Rapidity ordering:\n"; for(const auto & part: ev.outgoing()){ std::cout << std::setw(2) << part.type << ": "<< std::setw(7) << part.rapidity() << std::endl; } } int main(int argn, char** argv){ if(argn != 3){ std::cerr << "\n# Usage:\n."<< argv[0] <<" config.yml input_file.lhe\n\n"; return EXIT_FAILURE; } HEJ::Config config = HEJ::load_config(argv[1]); config.scales = HEJ::to_ScaleConfig( YAML::Load("scales: [H_T, 1 * H_T, 2/2 * H_T, 2*H_T/2, H_T/2*2, H_T/2/2*4, H_T*H_T/H_T]") ); HEJ::istream in{argv[2]}; LHEF::Reader reader{in}; auto ran = HEJ::make_RNG(config.rng.name, config.rng.seed); HEJ::ScaleGenerator scale_gen{ config.scales.base, config.scales.factors, config.scales.max_ratio }; HEJ::EventReweighter resum{ reader.heprup, std::move(scale_gen), to_EventReweighterConfig(config), *ran }; size_t i = 0; while(reader.readEvent()){ ++i; HEJ::Event event{ - HEJ::UnclusteredEvent{reader.hepeup}, + HEJ::EventData{reader.hepeup}, config.resummation_jets.def, config.resummation_jets.min_pt }; auto resummed = resum.reweight(event, config.trials); for(auto && ev: resummed) { for(auto &&var: ev.variations()) { if(std::abs(var.muf - ev.central().muf) > ep) { std::cerr << std::setprecision(15) << "unequal scales: " << var.muf << " != " << ev.central().muf << '\n' << "in resummed event:\n"; dump(ev); std::cerr << "\noriginal event:\n"; dump(event); return EXIT_FAILURE; } } } } } diff --git a/t/test_scale_import.cc b/t/test_scale_import.cc index 7a194c3..24ac7f3 100644 --- a/t/test_scale_import.cc +++ b/t/test_scale_import.cc @@ -1,31 +1,31 @@ #include #include #include "HEJ/YAMLreader.hh" #include "HEJ/Event.hh" int main(int argc, char** argv) { constexpr double ep = 1e-7; if (argc != 2) { throw std::logic_error{"wrong number of args"}; } const HEJ::Config config = HEJ::load_config(argv[1]); - HEJ::UnclusteredEvent tmp; - tmp.outgoing.push_back( + HEJ::EventData tmp; + tmp.outgoing().push_back( {HEJ::ParticleID::gluon, fastjet::PtYPhiM(50., -1., 0.3, 0.)} ); - tmp.outgoing.push_back( + tmp.outgoing().push_back( {HEJ::ParticleID::gluon, fastjet::PtYPhiM(30., 1., -0.3, 0.)} ); HEJ::Event ev{ std::move(tmp), fastjet::JetDefinition{fastjet::kt_algorithm, 0.4}, 20. }; const double softest_pt = config.scales.base[0](ev); if(std::abs(softest_pt-30.) > ep){ throw std::logic_error{"wrong softest pt"}; } }