diff --git a/include/HEJ/Config.hh b/include/HEJ/Config.hh index 12ac088..c134e23 100644 --- a/include/HEJ/Config.hh +++ b/include/HEJ/Config.hh @@ -1,234 +1,236 @@ /** \file * \brief HEJ 2 configuration parameters * * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #pragma once #include #include #include #include #include "fastjet/JetDefinition.hh" #include "yaml-cpp/yaml.h" #include "HEJ/Constants.hh" #include "HEJ/EWConstants.hh" #include "HEJ/Fraction.hh" #include "HEJ/HiggsCouplingSettings.hh" #include "HEJ/ScaleFunction.hh" #include "HEJ/event_types.hh" #include "HEJ/optional.hh" #include "HEJ/output_formats.hh" namespace HEJ { //! Jet parameters struct JetParameters{ fastjet::JetDefinition def; /**< Jet Definition */ double min_pt{}; /**< Minimum Jet Transverse Momentum */ }; //! Settings for scale variation struct ScaleConfig{ //! Base scale choices std::vector base; //! Factors for multiplicative scale variation std::vector factors; //! Maximum ratio between renormalisation and factorisation scale double max_ratio{}; }; //! Settings for random number generator struct RNGConfig { //! Random number generator name std::string name; //! Optional initial seed optional seed; }; //! Settings for partial unweighting struct PartialUnweightConfig { //! Number of trials for training size_t trials; //! Maximum distance in standard deviations from mean logarithmic weight double max_dev; }; /**! Possible treatments for fixed-order input events. * * The program will decide on how to treat an event based on * the value of this enumeration. */ enum class EventTreatment{ reweight, /**< Perform resummation */ keep, /**< Keep the event */ discard, /**< Discard the event */ }; //! Container to store the treatments for various event types using EventTreatMap = std::map; //! Possible setting for the event weight enum class WeightType{ weighted, //!< weighted events unweighted_resum, //!< unweighted only resummation part partially_unweighted //!< mixed weighted and unweighted }; /**! Input parameters. * * This struct handles stores all configuration parameters * needed in a HEJ 2 run. * * \internal To add a new option: * 1. Add a member to the Config struct. * 2. Inside "src/YAMLreader.cc": * - Add the option name to the "supported" Node in * get_supported_options. * - Initialise the new Config member in to_Config. * The functions set_from_yaml (for mandatory options) and * set_from_yaml_if_defined (non-mandatory) may be helpful. * 3. Add a new entry (with short description) to config.yaml * 4. Update the user documentation in "doc/Sphinx/" */ struct Config { //! %Parameters for scale variation ScaleConfig scales; //! Resummation jet properties JetParameters resummation_jets; //! Fixed-order jet properties JetParameters fixed_order_jets; //! Minimum transverse momentum for extremal partons //! \deprecated This will be removed in future versions. - //! Use \ref max_ext_soft_pt_fraction instead. - double min_extparton_pt = 0; - //! Maximum transverse momentum fraction from soft radiation in extremal jets - Fraction max_ext_soft_pt_fraction{}; + //! Use \ref max_regularisation_pt_fraction instead. + double min_extparton_pt = 0.; + //! \deprecated This is equivalent to\ref max_regularisation_pt_fraction + //! and will be removed in future versions. + optional> max_ext_soft_pt_fraction{}; //! The regulator lambda for the subtraction terms double regulator_lambda = CLAMBDA; //! Number of resummation configurations to generate per fixed-order event size_t trials{}; //! Maximal number of events optional max_events; //! Whether to include the logarithmic correction from \f$\alpha_s\f$ running bool log_correction{}; //! Event output files names and formats std::vector output; //! Parameters for random number generation RNGConfig rng; //! Map to decide what to do for different event types EventTreatMap treat; //! %Parameters for custom analysis //! @deprecated use analyses_parameters instead YAML::Node analysis_parameters; //! %Parameters for custom analyses std::vector analyses_parameters; //! Settings for effective Higgs-gluon coupling HiggsCouplingSettings Higgs_coupling; //! elector weak parameters EWConstants ew_parameters; //! Type of event weight e.g. (un)weighted WeightType weight_type; //! Settings for partial unweighting optional unweight_config; }; //! Configuration options for the PhaseSpacePoint class struct PhaseSpacePointConfig { //! Properties of resummation jets JetParameters jet_param; //! Minimum transverse momentum for extremal partons //! \deprecated This will be removed in future versions. - //! Use \ref max_ext_soft_pt_fraction instead. + //! Use \ref max_regularisation_pt_fraction instead. double min_extparton_pt = 0; - //! Maximum transverse momentum fraction from soft radiation in extremal jets - Fraction max_ext_soft_pt_fraction; + //! \deprecated This is equivalent to\ref max_regularisation_pt_fraction + //! and will be removed in future versions. + optional> max_ext_soft_pt_fraction{}; }; //! Configuration options for the MatrixElement class struct MatrixElementConfig { MatrixElementConfig() = default; MatrixElementConfig( bool log_correction, HiggsCouplingSettings Higgs_coupling, EWConstants ew_parameters, double regulator_lambda = CLAMBDA ): log_correction{log_correction}, Higgs_coupling{std::move(Higgs_coupling)}, ew_parameters{std::move(ew_parameters)}, regulator_lambda{regulator_lambda} {} //! Whether to include the logarithmic correction from \f$\alpha_s\f$ running bool log_correction{}; //! Settings for effective Higgs-gluon coupling HiggsCouplingSettings Higgs_coupling; //! elector weak parameters EWConstants ew_parameters; //! The regulator lambda for the subtraction terms double regulator_lambda = CLAMBDA; }; //! Configuration options for the EventReweighter class struct EventReweighterConfig { //! Settings for phase space point generation PhaseSpacePointConfig psp_config; //! Settings for matrix element calculation MatrixElementConfig ME_config; //! Access properties of resummation jets JetParameters & jet_param() { return psp_config.jet_param;} //! Access properties of resummation jets (const version) JetParameters const & jet_param() const { return psp_config.jet_param;} //! Treatment of the various event types EventTreatMap treat; }; /**! Extract PhaseSpacePointConfig from Config * * \internal We do not provide a PhaseSpacePointConfig constructor from Config * so that PhaseSpacePointConfig remains an aggregate. * This faciliates writing client code (e.g. the HEJ fixed-order generator) * that creates a PhaseSpacePointConfig *without* a Config object. * * @see to_MatrixElementConfig, to_EventReweighterConfig */ inline PhaseSpacePointConfig to_PhaseSpacePointConfig(Config const & conf) { return { conf.resummation_jets, conf.min_extparton_pt, conf.max_ext_soft_pt_fraction }; } /**! Extract MatrixElementConfig from Config * * @see to_PhaseSpacePointConfig, to_EventReweighterConfig */ inline MatrixElementConfig to_MatrixElementConfig(Config const & conf) { return {conf.log_correction, conf.Higgs_coupling, conf.ew_parameters, conf.regulator_lambda}; } /**! Extract EventReweighterConfig from Config * * @see to_PhaseSpacePointConfig, to_MatrixElementConfig */ inline EventReweighterConfig to_EventReweighterConfig(Config const & conf) { return { to_PhaseSpacePointConfig(conf), to_MatrixElementConfig(conf), conf.treat }; } } // namespace HEJ diff --git a/src/EventReweighter.cc b/src/EventReweighter.cc index b09dabf..17b725a 100644 --- a/src/EventReweighter.cc +++ b/src/EventReweighter.cc @@ -1,257 +1,257 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "HEJ/EventReweighter.hh" #include #include #include #include #include #include #include #include #include "fastjet/ClusterSequence.hh" #include "fastjet/PseudoJet.hh" #include "LHEF/LHEF.h" #include "HEJ/Event.hh" #include "HEJ/Fraction.hh" #include "HEJ/PDG_codes.hh" #include "HEJ/Particle.hh" #include "HEJ/PhaseSpacePoint.hh" #include "HEJ/exceptions.hh" namespace HEJ { EventReweighter::EventReweighter( LHEF::HEPRUP const & heprup, ScaleGenerator scale_gen, EventReweighterConfig conf, std::shared_ptr ran ): EventReweighter{ Beam{ heprup.EBMUP.first, {{ static_cast(heprup.IDBMUP.first), static_cast(heprup.IDBMUP.second) }} }, heprup.PDFSUP.first, std::move(scale_gen), std::move(conf), std::move(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 const & beam, int pdf_id, ScaleGenerator scale_gen, EventReweighterConfig conf, std::shared_ptr 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_{std::move(ran)} { assert(ran_); } PDF const & EventReweighter::pdf() const{ return pdf_; } std::vector EventReweighter::reweight( Event const & input_ev, std::size_t 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_(std::move(event)); return rescale(input_ev, std::move(res_events)); } EventTreatment EventReweighter::treatment(EventType type) const { return param_.treat.at(type); } std::vector EventReweighter::gen_res_events( Event const & ev, std::size_t phase_space_points ){ assert(ev.variations().empty()); status_.clear(); switch(treatment(ev.type())){ case EventTreatment::discard: { status_.emplace_back(StatusCode::discard); return {}; } case EventTreatment::keep: if(! jets_pass_resummation_cuts(ev)) { status_.emplace_back(StatusCode::failed_resummation_cuts); return {}; } else { status_.emplace_back(StatusCode::good); return {ev}; } default:; } const double Born_shat = shat(ev); std::vector resummation_events; status_.reserve(phase_space_points); for(std::size_t psp_number = 0; psp_number < phase_space_points; ++psp_number){ PhaseSpacePoint psp{ev, param_.psp_config, *ran_}; status_.emplace_back(psp.status()); assert(psp.status() != StatusCode::unspecified); if(psp.status() != StatusCode::good) continue; assert(psp.weight() != 0.); if(psp.incoming()[0].E() > E_beam_ || psp.incoming()[1].E() > E_beam_) { status_.back() = StatusCode::too_much_energy; continue; } resummation_events.emplace_back( to_EventData( std::move(psp) ).cluster( param_.jet_param().def, param_.jet_param().min_pt ) ); auto & new_event = resummation_events.back(); assert( new_event.valid_hej_state( - param_.psp_config.max_ext_soft_pt_fraction, + *param_.psp_config.max_ext_soft_pt_fraction, param_.psp_config.min_extparton_pt ) ); if( new_event.type() != ev.type() ) throw std::logic_error{"Resummation Event does not match Born event"}; new_event.generate_colours(*ran_); 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.parameters() *= pdf*ME/(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 if(part.type == pid::Higgs) { alpha_s_power += 2; } // nothing to do for other uncoloured particles } Weights result; result.central = std::pow(pdf_.Halphas(ev.central().mur), alpha_s_power); for(auto const & var: ev.variations()){ result.variations.emplace_back( std::pow(pdf_.Halphas(var.mur), alpha_s_power) ); } return result; } } // namespace HEJ diff --git a/src/PhaseSpacePoint.cc b/src/PhaseSpacePoint.cc index 81bd814..37415e6 100644 --- a/src/PhaseSpacePoint.cc +++ b/src/PhaseSpacePoint.cc @@ -1,856 +1,858 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "HEJ/PhaseSpacePoint.hh" #include #include #include #include #include #include #include #include #include #include #include #include "fastjet/ClusterSequence.hh" #include "fastjet/JetDefinition.hh" #include "HEJ/Constants.hh" #include "HEJ/Event.hh" #include "HEJ/JetSplitter.hh" #include "HEJ/PDG_codes.hh" #include "HEJ/RNG.hh" #include "HEJ/event_types.hh" #include "HEJ/kinematics.hh" #include "HEJ/resummation_jet.hh" #include "HEJ/utility.hh" namespace HEJ { 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; } namespace user_idx { //! user indices for partons with extremal rapidity enum ID: int { qqxmid1 = -9, qqxmid2 = -8, qqxb = -7, qqxf = -6, unob = -5, unof = -4, backward_fkl = -3, forward_fkl = -2, }; } // namespace user_idx using UID = user_idx::ID; double estimate_ng_mean(std::vector const & Born_jets){ const double delta_y = Born_jets.back().rapidity() - Born_jets.front().rapidity(); assert(delta_y > 0); // Formula derived from fit in arXiv:1805.04446 (see Fig. 2) return 0.975052*delta_y; } double phase_space_normalisation( int num_Born_jets, int num_out_partons ){ return std::pow(16.*std::pow(M_PI,3), num_Born_jets - num_out_partons); } } // namespace Event::EventData to_EventData(PhaseSpacePoint psp){ Event::EventData result; result.incoming = std::move(psp).incoming_; // NOLINT(bugprone-use-after-move) result.outgoing = std::move(psp).outgoing_; // NOLINT(bugprone-use-after-move) // technically Event::EventData doesn't have to be sorted, // but PhaseSpacePoint should be anyway assert( std::is_sorted( begin(result.outgoing), end(result.outgoing), rapidity_less{} ) ); assert(result.outgoing.size() >= 2); static_assert( std::numeric_limits::has_quiet_NaN, "no quiet NaN for double" ); constexpr double nan = std::numeric_limits::quiet_NaN(); result.decays = std::move(psp).decays_; // NOLINT(bugprone-use-after-move) result.parameters.central = {nan, nan, psp.weight()}; // NOLINT(bugprone-use-after-move) return result; } std::vector PhaseSpacePoint::cluster_jets( std::vector const & partons ) const{ fastjet::ClusterSequence cs(partons, param_.jet_param.def); return sorted_by_rapidity(cs.inclusive_jets(param_.jet_param.min_pt)); } bool PhaseSpacePoint::pass_resummation_cuts( std::vector const & jets ) const{ return cluster_jets(jets).size() == jets.size(); } int PhaseSpacePoint::sample_ng( std::vector const & Born_jets, RNG & ran ){ const double ng_mean = estimate_ng_mean(Born_jets); std::poisson_distribution dist(ng_mean); const int ng = dist(ran); assert(ng >= 0); assert(ng < NG_MAX); weight_ *= std::tgamma(ng + 1)*std::exp(ng_mean)*std::pow(ng_mean, -ng); return ng; } void PhaseSpacePoint::copy_AWZH_boson_from(Event const & event){ auto const & from = event.outgoing(); const auto AWZH_boson = std::find_if( begin(from), end(from), [](Particle const & p){ return is_AWZH_boson(p); } ); if(AWZH_boson == end(from)) return; auto insertion_point = std::lower_bound( begin(outgoing_), end(outgoing_), *AWZH_boson, rapidity_less{} ); outgoing_.insert(insertion_point, *AWZH_boson); // copy decay products const int idx = std::distance(begin(from), AWZH_boson); assert(idx >= 0); const auto decay_it = event.decays().find(idx); if(decay_it != end(event.decays())){ const int new_idx = std::distance(begin(outgoing_), insertion_point); assert(new_idx >= 0); assert(outgoing_[new_idx].type == AWZH_boson->type); decays_.emplace(new_idx, decay_it->second); } assert(std::is_sorted(begin(outgoing_), end(outgoing_), rapidity_less{})); } namespace { auto get_first_anyquark_emission(Event const & ev) { // find born quarks (ignore extremal partons) auto const firstquark = std::find_if( std::next(ev.begin_partons()), std::prev(ev.end_partons(), 2), [](Particle const & s){ return (is_anyquark(s)); } ); // assert that it is a q-q_bar pair. assert(std::distance(firstquark, ev.end_partons()) != 2); assert( ( is_quark(*firstquark) && is_antiquark(*std::next(firstquark)) ) || ( is_antiquark(*firstquark) && is_quark(*std::next(firstquark)) ) ); return firstquark; } //! returns index of most backward q-qbar jet template int get_back_quark_jet(Event const & ev, Iterator firstquark){ // find jets at FO corresponding to the quarks // technically this isn't necessary for LO std::vector const born_indices{ ev.particle_jet_indices() }; const auto firstquark_idx = std::distance(ev.begin_partons(), firstquark); int const firstjet_idx = born_indices[firstquark_idx]; assert(firstjet_idx>0); assert( born_indices[firstquark_idx+1] == firstjet_idx+1 ); return firstjet_idx; } //! returns index of most backward q-qbar jet int getBackQuarkJet(Event const & ev){ const auto firstquark = get_first_anyquark_emission(ev); return get_back_quark_jet(ev, firstquark); } template void label_extremal_qqx( ConstIterator born_begin, ConstIterator born_end, Iterator first_out ){ // find born quarks const auto firstquark = std::find_if( born_begin, born_end-1, [](Particle const & s){ return (is_anyquark(s)); } ); assert(firstquark != born_end-1); const auto secondquark = std::find_if( firstquark+1, born_end, [](Particle const & s){ return (is_anyquark(s)); } ); assert(secondquark != born_end); assert( ( is_quark(*firstquark) && is_antiquark(*secondquark) ) || ( is_antiquark(*firstquark) && is_quark(*secondquark) )); assert(first_out->type == ParticleID::gluon); assert((first_out+1)->type == ParticleID::gluon); // copy type from born first_out->type = firstquark->type; (first_out+1)->type = secondquark->type; } } // namespace void PhaseSpacePoint::label_qqx(Event const & event){ assert(std::is_sorted(begin(outgoing_), end(outgoing_), rapidity_less{})); assert(filter_partons(outgoing_).size() == outgoing_.size()); if(qqxb_){ label_extremal_qqx( event.outgoing().cbegin(), event.outgoing().cend(), outgoing_.begin() ); return; } if(qqxf_){ // same as qqxb with reversed order label_extremal_qqx( event.outgoing().crbegin(), event.outgoing().crend(), outgoing_.rbegin() ); return; } // central qqx const auto firstquark = get_first_anyquark_emission(event); // find jets at FO corresponding to the quarks // technically this isn't necessary for LO const auto firstjet_idx = get_back_quark_jet(event, firstquark); // find corresponding jets after resummation fastjet::ClusterSequence cs{to_PseudoJet(outgoing_), param_.jet_param.def}; auto const jets = fastjet::sorted_by_rapidity( cs.inclusive_jets( param_.jet_param.min_pt )); std::vector const resum_indices{ cs.particle_jet_indices({jets}) }; // assert that jets didn't move assert(nearby_ep( ( event.jets().cbegin()+firstjet_idx )->rapidity(), jets[ firstjet_idx ].rapidity(), 1e-2) ); assert(nearby_ep( ( event.jets().cbegin()+firstjet_idx+1 )->rapidity(), jets[ firstjet_idx+1 ].rapidity(), 1e-2) ); // find last partons in first (central) jet size_t idx_out = 0; for(size_t i=resum_indices.size()-2; i>0; --i) if(resum_indices[i] == firstjet_idx){ idx_out = i; break; } assert(idx_out != 0); // check that there is sufficient pt in jets from the quarks - const double minpartonjetpt = 1. - param_.max_ext_soft_pt_fraction; + const double minpartonjetpt = 1. - *param_.max_ext_soft_pt_fraction; if (outgoing_[idx_out].p.pt()pt()){ weight_=0.; status_ = StatusCode::wrong_jets; return; } if (outgoing_[idx_out+1].p.pt()pt()){ weight_=0.; status_ = StatusCode::wrong_jets; return; } // check that no additional emission between jets // such configurations are possible if we have an gluon gets generated // inside the rapidities of the qqx chain, but clusted to a // differnet/outside jet. Changing this is non trivial if(resum_indices[idx_out+1] != resum_indices[idx_out]+1){ weight_=0.; status_ = StatusCode::gluon_in_qqx; return; } outgoing_[idx_out].type = firstquark->type; outgoing_[idx_out+1].type = std::next(firstquark)->type; } void PhaseSpacePoint::label_quarks(Event const & ev){ const auto WZEmit = std::find_if( begin(ev.outgoing()), end(ev.outgoing()), [](Particle const & s){ return (std::abs(s.type) == pid::Wp || s.type == pid::Z_photon_mix); } ); if (WZEmit != end(ev.outgoing())){ if(!qqxb_) { const size_t backward_FKL_idx = unob_?1:0; const auto backward_FKL = std::next(ev.begin_partons(), backward_FKL_idx); outgoing_[backward_FKL_idx].type = backward_FKL->type; } if(!qqxf_) { const size_t forward_FKL_idx = unof_?1:0; const auto forward_FKL = std::prev(ev.end_partons(), 1+forward_FKL_idx); outgoing_.rbegin()[unof_].type = forward_FKL->type; // NOLINT } } else { most_backward_FKL(outgoing_).type = ev.incoming().front().type; most_forward_FKL(outgoing_).type = ev.incoming().back().type; } if(qqxmid_||qqxb_||qqxf_){ label_qqx(ev); } } PhaseSpacePoint::PhaseSpacePoint( Event const & ev, PhaseSpacePointConfig conf, RNG & ran ): unob_{ev.type() == event_type::unob}, unof_{ev.type() == event_type::unof}, qqxb_{ev.type() == event_type::qqxexb}, qqxf_{ev.type() == event_type::qqxexf}, qqxmid_{ev.type() == event_type::qqxmid}, param_{std::move(conf)}, status_{unspecified} { + if(!param_.max_ext_soft_pt_fraction) + throw HEJ::missing_option{"\"max_ext_soft_pt_fraction\" not set!"}; weight_ = 1; auto const & Born_jets = ev.jets(); const int ng = sample_ng(Born_jets, ran); weight_ /= std::tgamma(ng + 1); const int ng_jets = sample_ng_jets(ng, Born_jets, ran); std::vector out_partons = gen_non_jet( ng - ng_jets, CMINPT, param_.jet_param.min_pt, ran ); int qqxbackjet(-1); if(qqxmid_){ qqxbackjet = getBackQuarkJet(ev); } const auto qperp = std::accumulate( begin(out_partons), end(out_partons), fastjet::PseudoJet{} ); const auto jets = reshuffle(Born_jets, qperp); if(weight_ == 0.) { status_ = failed_reshuffle; return; } if(! pass_resummation_cuts(jets)){ status_ = failed_resummation_cuts; weight_ = 0.; return; } std::vector jet_partons = split( jets, ng_jets, qqxbackjet, ran ); if(weight_ == 0.) { status_ = StatusCode::failed_split; return; } if(qqxmid_){ rescale_qqx_rapidities( out_partons, jets, most_backward_FKL(jet_partons).rapidity(), most_forward_FKL(jet_partons).rapidity(), qqxbackjet ); } else{ rescale_rapidities( out_partons, most_backward_FKL(jet_partons).rapidity(), most_forward_FKL(jet_partons).rapidity() ); } if(! cluster_jets(out_partons).empty()){ weight_ = 0.; status_ = StatusCode::empty_jets; return; } std::sort(begin(out_partons), end(out_partons), rapidity_less{}); assert( std::is_sorted(begin(jet_partons), end(jet_partons), rapidity_less{}) ); const auto first_jet_parton = out_partons.insert( end(out_partons), begin(jet_partons), end(jet_partons) ); std::inplace_merge( begin(out_partons), first_jet_parton, end(out_partons), rapidity_less{} ); if(! jets_ok(Born_jets, out_partons)){ weight_ = 0.; status_ = StatusCode::wrong_jets; return; } weight_ *= phase_space_normalisation(Born_jets.size(), out_partons.size()); outgoing_.reserve(out_partons.size() + 1); // one slot for possible A, W, Z, H for( auto it = std::make_move_iterator(out_partons.begin()); it != std::make_move_iterator(out_partons.end()); ++it ){ outgoing_.emplace_back( Particle{pid::gluon, *it, {}}); } assert(!outgoing_.empty()); label_quarks(ev); if(weight_ == 0.) { //! @TODO optimise s.t. this is not possible // status is handled internally return; } copy_AWZH_boson_from(ev); reconstruct_incoming(ev.incoming()); status_ = StatusCode::good; } std::vector PhaseSpacePoint::gen_non_jet( int const ng_non_jet, double const ptmin, double const ptmax, RNG & ran ){ // heuristic parameters for pt sampling const double ptpar = 1.3 + ng_non_jet/5.; const double temp1 = std::atan((ptmax - ptmin)/ptpar); std::vector partons(ng_non_jet); for(int i = 0; i < ng_non_jet; ++i){ const double r1 = ran.flat(); const double pt = ptmin + ptpar*std::tan(r1*temp1); const double temp2 = std::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 = ran.flat(); 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 + NG_MAX); assert(ptmin-1e-5 <= partons[i].pt() && partons[i].pt() <= ptmax+1e-5); } assert(std::all_of(partons.cbegin(), partons.cend(), is_nonjet_parton)); return sorted_by_rapidity(partons); } void PhaseSpacePoint::rescale_qqx_rapidities( std::vector & out_partons, std::vector const & jets, const double ymin1, const double ymax2, const int qqxbackjet ){ const double ymax1 = jets[qqxbackjet].rapidity(); const double ymin2 = jets[qqxbackjet+1].rapidity(); constexpr double ep = 1e-7; const double tot_y = ymax1 - ymin1 + ymax2 - ymin2; std::vector> refpart( out_partons.begin(), out_partons.end()); double ratio = (ymax1 - ymin1)/tot_y; const auto gap{ std::find_if(refpart.begin(), refpart.end(), [ratio](fastjet::PseudoJet const & p){ return (p.rapidity()>=ratio);} ) }; double ymin = ymin1; double ymax = ymax1; double dy = ymax - ymin - 2*ep; double offset = 0.; for(auto it_part=refpart.begin(); it_part & 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); } } namespace { template auto min(T const & a, T const & b, Rest&&... r) { using std::min; return min(a, min(b, std::forward(r)...)); } } double PhaseSpacePoint::probability_in_jet( std::vector const & Born_jets ) const{ assert(std::is_sorted(begin(Born_jets), end(Born_jets), rapidity_less{})); assert(Born_jets.size() >= 2); const double dy = Born_jets.back().rapidity() - Born_jets.front().rapidity(); const double R = param_.jet_param.def.R(); const int njets = Born_jets.size(); const double p_J_y_large = (njets-1)*R*R/(2.*dy); const double p_J_y0 = njets*R/M_PI; return min(p_J_y_large, p_J_y0, 1.); } int PhaseSpacePoint::sample_ng_jets( int ng, std::vector const & Born_jets, 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(ran); weight_ *= std::pow(p_J, -ng_J)*std::pow(1 - p_J, ng_J - ng); return ng_J; } std::vector PhaseSpacePoint::reshuffle( std::vector const & Born_jets, fastjet::PseudoJet const & q ){ if(q == fastjet::PseudoJet{0, 0, 0, 0}) return Born_jets; auto jets = resummation_jet_momenta(Born_jets, q); if(jets.empty()){ weight_ = 0; return {}; } // additional Jacobian to ensure Born integration over delta gives 1 weight_ *= resummation_jet_weight(Born_jets, q); return jets; } std::vector PhaseSpacePoint::distribute_jet_partons( int ng_jets, std::vector const & jets, RNG & ran ){ size_t first_valid_jet = 0; size_t num_valid_jets = jets.size(); const double R_eff = 5./3.*param_.jet_param.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_||qqxb_) && jets[0].delta_R(jets[1]) > R_eff){ ++first_valid_jet; --num_valid_jets; } else if((unof_||qqxf_) && jets[jets.size()-1].delta_R(jets[jets.size()-2]) > R_eff){ --num_valid_jets; } std::vector 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; } #ifndef NDEBUG namespace { bool tagged_FKL_backward( std::vector const & jet_partons ){ return std::find_if( begin(jet_partons), end(jet_partons), [](fastjet::PseudoJet const & p){ return p.user_index() == UID::backward_fkl; } ) != end(jet_partons); } bool tagged_FKL_forward( std::vector const & jet_partons ){ // the most forward FKL parton is most likely near the end of jet_partons; // start search from there return std::find_if( jet_partons.rbegin(), jet_partons.rend(), [](fastjet::PseudoJet const & p){ return p.user_index() == UID::forward_fkl; } ) != jet_partons.rend(); } bool tagged_FKL_extremal( std::vector const & jet_partons ){ return tagged_FKL_backward(jet_partons) && tagged_FKL_forward(jet_partons); } } // namespace #endif std::vector PhaseSpacePoint::split( std::vector const & jets, int ng_jets, size_t qqxbackjet, RNG & ran ){ return split( jets, distribute_jet_partons(ng_jets, jets, ran), qqxbackjet, ran); } bool PhaseSpacePoint::pass_extremal_cuts( fastjet::PseudoJet const & ext_parton, fastjet::PseudoJet const & jet ) const{ if(ext_parton.pt() < param_.min_extparton_pt) return false; - return (ext_parton - jet).pt()/jet.pt() < param_.max_ext_soft_pt_fraction; + return (ext_parton - jet).pt()/jet.pt() < *param_.max_ext_soft_pt_fraction; } std::vector PhaseSpacePoint::split( std::vector const & jets, std::vector const & np, size_t qqxbackjet, RNG & ran ){ assert(! jets.empty()); assert(jets.size() == np.size()); assert(pass_resummation_cuts(jets)); const size_t most_backward_FKL_idx = 0 + unob_ + qqxb_; // NOLINT const size_t most_forward_FKL_idx = jets.size() - 1 - unof_ - qqxf_; // NOLINT auto const & jet = param_.jet_param; const JetSplitter jet_splitter{jet.def, jet.min_pt}; std::vector jet_partons; // randomly distribute jet gluons among jets for(size_t i = 0; i < jets.size(); ++i){ auto split_res = jet_splitter.split(jets[i], np[i], ran); weight_ *= split_res.weight; if(weight_ == 0) return {}; assert( std::all_of( begin(split_res.constituents), end(split_res.constituents), is_jet_parton ) ); const auto first_new_parton = jet_partons.insert( end(jet_partons), begin(split_res.constituents), end(split_res.constituents) ); // mark uno and extremal FKL emissions here so we can check // their position once all emissions are generated // also mark qqxmid partons, and apply appropriate pt cut. auto extremal = end(jet_partons); if (i == most_backward_FKL_idx){ //FKL backward emission extremal = std::min_element( first_new_parton, end(jet_partons), rapidity_less{} ); extremal->set_user_index(UID::backward_fkl); } else if(((unob_ || qqxb_) && i == 0)){ // unordered/qqxb extremal = std::min_element( first_new_parton, end(jet_partons), rapidity_less{} ); extremal->set_user_index((unob_)?UID::unob:UID::qqxb); } else if (i == most_forward_FKL_idx){ extremal = std::max_element( first_new_parton, end(jet_partons), rapidity_less{} ); extremal->set_user_index(UID::forward_fkl); } else if(((unof_ || qqxf_) && i == jets.size() - 1)){ // unordered/qqxf extremal = std::max_element( first_new_parton, end(jet_partons), rapidity_less{} ); extremal->set_user_index((unof_)?UID::unof:UID::qqxf); } else if((qqxmid_ && i == qqxbackjet)){ extremal = std::max_element( first_new_parton, end(jet_partons), rapidity_less{} ); extremal->set_user_index(UID::qqxmid1); } else if((qqxmid_ && i == qqxbackjet+1)){ extremal = std::min_element( first_new_parton, end(jet_partons), rapidity_less{} ); extremal->set_user_index(UID::qqxmid2); } 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_ok(jet_partons) || !split_preserved_jets(jets, jet_partons) ){ weight_ = 0.; return {}; } return jet_partons; } bool PhaseSpacePoint::extremal_ok( std::vector const & partons ) const{ assert(std::is_sorted(begin(partons), end(partons), rapidity_less{})); if(unob_ && partons.front().user_index() != UID::unob) return false; if(unof_ && partons.back().user_index() != UID::unof) return false; if(qqxb_ && partons.front().user_index() != UID::qqxb) return false; if(qqxf_ && partons.back().user_index() != UID::qqxf) return false; return most_backward_FKL(partons).user_index() == UID::backward_fkl && most_forward_FKL(partons).user_index() == UID::forward_fkl; } bool PhaseSpacePoint::split_preserved_jets( std::vector const & jets, std::vector const & jet_partons ) const{ assert(std::is_sorted(begin(jets), end(jets), rapidity_less{})); const auto split_jets = 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 Particle const & PhaseSpacePoint::most_backward_FKL( std::vector const & partons ) const{ return partons[0 + unob_ + qqxb_]; } template Particle const & PhaseSpacePoint::most_forward_FKL( std::vector const & partons ) const{ const size_t idx = partons.size() - 1 - unof_ - qqxf_; assert(idx < partons.size()); return partons[idx]; } template Particle & PhaseSpacePoint::most_backward_FKL( std::vector & partons ) const{ return partons[0 + unob_ + qqxb_]; } template Particle & PhaseSpacePoint::most_forward_FKL( std::vector & partons ) const{ const size_t idx = partons.size() - 1 - unof_ - qqxf_; assert(idx < partons.size()); return partons[idx]; } bool PhaseSpacePoint::contains_idx( fastjet::PseudoJet const & jet, fastjet::PseudoJet const & parton ) const { auto const & constituents = jet.constituents(); const int idx = parton.user_index(); const bool injet = std::find_if( begin(constituents), end(constituents), [idx](fastjet::PseudoJet const & con){return con.user_index() == idx;} ) != end(constituents); - const double minpartonjetpt = 1. - param_.max_ext_soft_pt_fraction; + const double minpartonjetpt = 1. - *param_.max_ext_soft_pt_fraction; return ((parton.pt()>minpartonjetpt*jet.pt())&&injet); } bool PhaseSpacePoint::jets_ok( std::vector const & Born_jets, std::vector const & partons ) const{ fastjet::ClusterSequence cs(partons, param_.jet_param.def); const auto jets = sorted_by_rapidity(cs.inclusive_jets(param_.jet_param.min_pt)); if(jets.size() != Born_jets.size()) return false; int in_jet = 0; for(auto const & jet : jets){ assert(jet.has_constituents()); for(auto && parton: jet.constituents()){ if(is_nonjet_parton(parton)) return false; } in_jet += jet.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; if(unob_ && !contains_idx(jets.front(), partons.front())) return false; if(qqxb_ && !contains_idx(jets.front(), partons.front())) return false; if(unof_ && !contains_idx(jets.back(), partons.back())) return false; if(qqxf_ && !contains_idx(jets.back(), partons.back())) return false; #ifndef NDEBUG for(size_t i = 0; i < jets.size(); ++i){ assert(nearby_ep(jets[i].rapidity(), Born_jets[i].rapidity(), 1e-2)); } #endif return true; } void PhaseSpacePoint::reconstruct_incoming( std::array 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()); } bool PhaseSpacePoint::momentum_conserved() const{ fastjet::PseudoJet diff; for(auto const & in: incoming()) diff += in.p; const double norm = diff.E(); for(auto const & out: outgoing()) diff -= out.p; return nearby(diff, fastjet::PseudoJet{}, norm); } } //namespace HEJ diff --git a/t/check_res.cc b/t/check_res.cc index 101a4fa..177d5e2 100644 --- a/t/check_res.cc +++ b/t/check_res.cc @@ -1,166 +1,166 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "hej_test.hh" #include #include #include #include #include #include #include #include #include "HEJ/Config.hh" #include "HEJ/CrossSectionAccumulator.hh" #include "HEJ/Event.hh" #include "HEJ/event_types.hh" #include "HEJ/EventReweighter.hh" #include "HEJ/EWConstants.hh" #include "HEJ/Fraction.hh" #include "HEJ/HiggsCouplingSettings.hh" #include "HEJ/Mixmax.hh" #include "HEJ/Parameters.hh" #include "HEJ/ScaleFunction.hh" #include "HEJ/stream.hh" #include "fastjet/JetDefinition.hh" #include "LHEF/LHEF.h" namespace HEJ { struct RNG; } namespace { const fastjet::JetDefinition JET_DEF{fastjet::kt_algorithm, 0.4}; const fastjet::JetDefinition BORN_JET_DEF{JET_DEF}; constexpr double BORN_JETPTMIN = 30; - constexpr double MAX_EXT_SOFT_PT_FRACTION = 0.1; + const HEJ::Fraction MAX_EXT_SOFT_PT_FRACTION{0.1}; constexpr double JETPTMIN = 35; constexpr bool LOG_CORR = false; constexpr std::size_t NUM_TRIES = 100; constexpr HEJ::ParticleProperties WPROP{80.385, 2.085}; constexpr HEJ::ParticleProperties ZPROP{91.187, 2.495}; constexpr HEJ::ParticleProperties HPROP{125, 0.004165}; constexpr double VEV = 246.2196508; using EventTreatment = HEJ::EventTreatment; using namespace HEJ::event_type; HEJ::EventTreatMap TREAT{ {no_2_jets, EventTreatment::discard}, {bad_final_state, EventTreatment::discard}, {non_resummable, EventTreatment::discard}, {unof, EventTreatment::discard}, {unob, EventTreatment::discard}, {qqxexb, EventTreatment::discard}, {qqxexf, EventTreatment::discard}, {qqxmid, EventTreatment::discard}, {FKL, EventTreatment::reweight} }; bool correct_colour(HEJ::Event const & ev){ if(!HEJ::event_type::is_resummable(ev.type())) return true; return ev.is_leading_colour(); } } // namespace int main(int argn, char** argv) { if(argn == 5 && std::string(argv[4]) == "unof"){ --argn; TREAT[unof] = EventTreatment::reweight; TREAT[unob] = EventTreatment::discard; TREAT[FKL] = EventTreatment::discard; } if(argn == 5 && std::string(argv[4]) == "unob"){ --argn; TREAT[unof] = EventTreatment::discard; TREAT[unob] = EventTreatment::reweight; TREAT[FKL] = EventTreatment::discard; } else if(argn == 5 && std::string(argv[4]) == "splitf"){ --argn; TREAT[qqxexb] = EventTreatment::discard; TREAT[qqxexf] = EventTreatment::reweight; TREAT[FKL] = EventTreatment::discard; } else if(argn == 5 && std::string(argv[4]) == "splitb"){ --argn; TREAT[qqxexb] = EventTreatment::reweight; TREAT[qqxexf] = EventTreatment::discard; TREAT[FKL] = EventTreatment::discard; } else if(argn == 5 && std::string(argv[4]) == "qqxmid"){ --argn; TREAT[qqxmid] = EventTreatment::reweight; TREAT[FKL] = EventTreatment::discard; } if(argn != 4){ std::cerr << "Usage: check_res eventfile xsection tolerance [uno]"; return EXIT_FAILURE; } const double xsec_ref = std::stod(argv[2]); const double tolerance = std::stod(argv[3]); HEJ::istream in{argv[1]}; LHEF::Reader reader{in}; HEJ::PhaseSpacePointConfig psp_conf; psp_conf.jet_param = HEJ::JetParameters{JET_DEF, JETPTMIN}; psp_conf.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{}; ME_conf.ew_parameters.set_vevWZH(VEV, WPROP, ZPROP, HPROP); HEJ::EventReweighterConfig conf; conf.psp_config = std::move(psp_conf); conf.ME_config = std::move(ME_conf); conf.treat = TREAT; reader.readEvent(); const bool has_Higgs = std::find( begin(reader.hepeup.IDUP), end(reader.hepeup.IDUP), 25 ) != end(reader.hepeup.IDUP); const double mu = has_Higgs?125.:91.188; HEJ::ScaleGenerator scale_gen{ {{std::to_string(mu), HEJ::FixedScale{mu}}}, {}, 1. }; std::shared_ptr ran{std::make_shared()}; HEJ::EventReweighter hej{reader.heprup, std::move(scale_gen), conf, ran}; HEJ::CrossSectionAccumulator xs; do{ auto ev_data = HEJ::Event::EventData{reader.hepeup}; shuffle_particles(ev_data); ev_data.reconstruct_intermediate(); HEJ::Event ev{ ev_data.cluster( BORN_JET_DEF, BORN_JETPTMIN ) }; auto resummed_events = hej.reweight(ev, NUM_TRIES); for(auto const & res_ev: resummed_events) { ASSERT(correct_colour(res_ev)); ASSERT(std::isfinite(res_ev.central().weight)); // we fill the xs uncorrelated since we only want to test the uncertainty // of the resummation xs.fill(res_ev); } } while(reader.readEvent()); const double xsec = xs.total().value; const double xsec_err = std::sqrt(xs.total().error); const double significance = std::abs(xsec - xsec_ref) / std::sqrt( xsec_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; } return EXIT_SUCCESS; } diff --git a/t/test_psp.cc b/t/test_psp.cc index 5abc716..14d9f23 100644 --- a/t/test_psp.cc +++ b/t/test_psp.cc @@ -1,74 +1,74 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "hej_test.hh" #include #include #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" #include "LHEF/LHEF.h" #include "fastjet/JetDefinition.hh" namespace { constexpr int MAX_TRIALS = 100; - constexpr double MAX_EXT_SOFT_PT_FRACTION = 0.1; + const HEJ::Fraction MAX_EXT_SOFT_PT_FRACTION{0.1}; 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: " << argv[0] << " 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.max_ext_soft_pt_fraction = MAX_EXT_SOFT_PT_FRACTION; HEJ::Ranlux64 ran{}; while(reader.readEvent()){ HEJ::Event::EventData ev_data{ reader.hepeup }; shuffle_particles(ev_data); const HEJ::Event ev{ ev_data( JET_DEF, MIN_JET_PT ) }; for(int trial = 0; trial < MAX_TRIALS; ++trial){ HEJ::PhaseSpacePoint psp{ev, conf, ran}; if(psp.weight() != 0){ HEJ::Event::EventData tmp_ev; tmp_ev.incoming = psp.incoming(); tmp_ev.outgoing = psp.outgoing(); tmp_ev.parameters.central = {0,0,0}; shuffle_particles(tmp_ev); HEJ::Event out_ev{ tmp_ev(JET_DEF, MIN_JET_PT) }; if(out_ev.type() != ev.type()){ using HEJ::event_type::name; std::cerr << "Wrong class of phase space point:\n" "original event of class " << name(ev.type()) << ":\n"; writer.hepeup = reader.hepeup; writer.writeEvent(); std::cerr << "Phase space point of class " << name(out_ev.type()) << ":\n"; writer.hepeup = to_HEPEUP(out_ev, &writer.heprup); writer.writeEvent(); return EXIT_FAILURE; } } } } return EXIT_SUCCESS; }