diff --git a/Changes.md b/Changes.md index e6d2ade..8312713 100644 --- a/Changes.md +++ b/Changes.md @@ -1,140 +1,140 @@ # Changelog This is the log for changes to the HEJ program. Further changes to the HEJ API are documented in [`Changes-API.md`](Changes-API.md). If you are using HEJ as a library, please also read the changes there. ## Version 2.1 This release adds support for two new processes: * W boson with jets. * Jets with a charged lepton-antilepton pair via a virtual Z boson or photon. In addition, the complete set of first subleading processes (unordered gluon, central and extremal quark-antiquark pair) is implemented for pure jets and W + jets, see [arXiv:TODO](https://arxiv.org/abs/TODO). Unordered gluon emission is also supported for Higgs boson + jets and Z boson/photon + jets. This release include many changes to the code, which affect users of HEJ as a library (see [`Changes-API.md`](Changes-API.md)). ### 2.1.0 #### New Processes * Resummation for W bosons with jets - New subleading processes `extremal qqx` & `central qqx` for a quark and anti-quark in the final state, e.g. `g g => u d_bar Wm g` (the other subleading processes also work with W's) - `HEJFOG` can generate multiple jets together with a (off-shell) W bosons decaying into lepton & neutrino * Resummation for jets with a charged lepton-antilepton pair via a virtual Z boson or photon. Includes the `unordered` subleading process. * Resummation can now be performed on all subleading processes within pure jets also. This includes `unordered`, `extremal qqbar` and `central qqbar` processes. #### More Physics implementation * Partons now have a Colour charge - Colours are read from and written to LHE files - For reweighted events the colours are created according to leading colour in the FKL limit * Use relative fraction for soft transverse momentum in extremal jets (`max ext soft pt fraction`). This supersedes `min extparton pt`, which is now deprecated and will be removed in future versions. #### Updates to Runcard * Allow multiplication and division of multiple scale functions e.g. `H_T/2*m_j1j2` * Grouped `event treatment` for subleading channels together in runcard - Rename `non-HEJ` processes to `non-resummable` * Read electro-weak constants from input - new mandatory setting `vev` to change vacuum expectation value - new mandatory settings `particle properties` to specify mass & width of bosons - - FOG: decays are now specified in `decays` setting (previously under + - `HEJFOG`: decays are now specified in `decays` setting (previously under `particle properties`) * Allow loading multiple analyses with `analyses`. The old `analysis` (with "i") is marked deprecated. * Optional setting to specify maximal number of Fixed Order events (`max events`, default is all) * Allow changing the regulator lambda in input (`regulator parameter`, only for advanced users) #### Changes to Input/Output * Added support to read & write `hdf5` event files suggested in [arXiv:1905.05120](https://arxiv.org/abs/1905.05120) (needs [HighFive](https://github.com/BlueBrain/HighFive)) * Support input with average weight equal to the cross section (`IDWTUP=1 or 4`) * Support unmodified Les Houches Event Files written by Sherpa with `cross section = sum(weights)/sum(trials)` * Analyses now get general run information (`LHEF::HEPRUP`) in the constructor. **This might break previously written, external analyses!** - external analyses should now be created with `make_analysis(YAML::Node const & config, LHEF::HEPRUP const & heprup)` * Support `rivet` version 3 with both `HepMC` version 2 and 3 - Multiple weights with `rivet 3` will only create one `.yoda` file (instead of one per weight/scale) * Added option to unweight only resummation events (`unweight: {type: resummation}`) * Added option for partially unweighting resummation events, similar to the fixed-order generator. * Improved unweighting algorithm. * Follow HepMC convention for particle Status codes: incoming = 11, decaying = 2, outgoing = 1 (unchanged) #### Miscellaneous * Print cross sections at end of run * Added example analysis & scale to `examples/`. Everything in `examples/` will be build when the flag `-DBUILD_EXAMPLES=TRUE` is set in `cmake`. * Dropped support for HepMC 3.0.0, either HepMC version 2 or >3.1 is required - It is now possible to write out both HepMC 2 and HepMC 3 events at the same time * Require LHADPF version 6. Dropped support for all other versions. * Use `git-lfs` for raw data in test (`make test` now requires `git-lfs`) * Currents are now generated with [`FORM`](https://github.com/vermaseren/form) - `FORM` is included as a `git submodule`, use `git submodule update --init` to download `FORM` * Create [Sphinx](http://sphinx-doc.org/) and [Doxygen](http://doxygen.org/) documentation by `make sphinx` or `make doxygen` in your `build/` folder ## Version 2.0 First release of HEJ 2. Complete code rewrite compared to HEJ 1. Improved matching to Fixed Order ([arXiv:1805.04446](https://arxiv.org/abs/1805.04446)). Implemented processes: Higgs boson with jets (FKL and unordered gluon emission, with finite quark mass loop, [arXiv:1812.08072](https://arxiv.org/abs/1812.08072)), and pure jets (only FKL). See [arXiv:1902.08430](https://arxiv.org/abs/1902.08430) ### 2.0.6 * Fixed compiling rivet when YODA headers are _outside_ of rivet directory ### 2.0.5 * Fixed event classification for input not ordered in rapidity ### 2.0.4 * Fixed wrong path of `HEJ_INCLUDE_DIR` in `hej-config.cmake` ### 2.0.3 * Fixed parsing of (numerical factor) * (base scale) in configuration * Don't change scale names, but sanitise Rivet output file names instead ### 2.0.2 * Changed scale names to `"_over_"` and `"_times_"` for proper file names (was `"/"` and `"*"` before) ### 2.0.1 * Fixed name of fixed-order generator in error message. ### 2.0.0 * First release diff --git a/FixedOrderGen/configFO.yml b/FixedOrderGen/configFO.yml index 6fca80b..2ffcb83 100644 --- a/FixedOrderGen/configFO.yml +++ b/FixedOrderGen/configFO.yml @@ -1,94 +1,95 @@ ## Number of generated events events: 200 jets: min pt: 20 # minimal jet pt, should be slightly below analysis cut peak pt: 30 # peak pt of jets, should be at analysis cut algorithm: antikt # jet clustering algorithm R: 0.4 # jet R parameter max rapidity: 5 # maximum jet rapidity ## Particle beam beam: energy: 6500 particles: [p, p] ## PDF ID according to LHAPDF pdf: 230000 ## Scattering process process: p p => h 4j ## Particle decays (multiple decays are allowed) decays: Higgs: {into: [photon, photon], branching ratio: 0.0023568762400521404} Wp: {into: [e+, nu_e], branching ratio: 1} Wm: {into: [e-, nu_e_bar]} # equivalent to "branching ratio: 1" ## Fraction of events with two extremal emissions in one direction ## that contain an subleading emission e.g. unordered emission subleading fraction: 0.05 ## Allow different subleading configurations. ## By default all processes are allowed. ## This does not check if the processes are implemented in HEJ! # # subleading channels: # - unordered -# - qqx +# - central qqx +# - extremal qqx ## Unweighting parameters ## remove to obtain weighted events # unweight: # sample size: 200 # should be similar to "events:", but not more than ~10000 # max deviation: 0 ## --- The following settings are equivalent to HEJ, --- ## ## --- see HEJ for reference. --- ## ## Central scale choice scales: max jet pperp ## Event output files event output: - HEJFO.lhe # - HEJFO_events.hepmc ## Analyses # # analyses: ## Rivet analysis # - rivet: MC_XS # rivet analysis name # output: HEJ # name of the yoda files, ".yoda" and scale suffix will be added ## Custom analysis # - plugin: /path/to/libmyanalysis.so # my analysis parameter: some value ## Selection of random number generator and seed random generator: name: mixmax # seed: 1 ## Vacuum expectation value vev: 246.2196508 ## Properties of the weak gauge bosons particle properties: Higgs: mass: 125 width: 0.004165 W: mass: 80.385 width: 2.085 Z: mass: 91.187 width: 2.495 ## Parameters for Higgs-gluon couplings ## this requires compilation with QCDloop # # Higgs coupling: # use impact factors: false # mt: 174 # include bottom: true # mb: 4.7 diff --git a/FixedOrderGen/include/PhaseSpacePoint.hh b/FixedOrderGen/include/PhaseSpacePoint.hh index f35c164..d0224b1 100644 --- a/FixedOrderGen/include/PhaseSpacePoint.hh +++ b/FixedOrderGen/include/PhaseSpacePoint.hh @@ -1,282 +1,281 @@ /** \file PhaseSpacePoint.hh * \brief Contains the PhaseSpacePoint Class * * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #pragma once #include #include #include #include #include #include "boost/iterator/filter_iterator.hpp" #include "HEJ/Event.hh" #include "HEJ/Particle.hh" #include "HEJ/PDG_codes.hh" #include "fastjet/PseudoJet.hh" #include "Decay.hh" #include "Status.hh" #include "Subleading.hh" namespace HEJ { class EWConstants; class PDF; class RNG; } namespace HEJFOG { class JetParameters; class Process; //! A point in resummation phase space class PhaseSpacePoint{ public: //! Iterator over partons using ConstPartonIterator = boost::filter_iterator< bool (*)(HEJ::Particle const &), std::vector::const_iterator >; //! Reverse Iterator over partons using ConstReversePartonIterator = std::reverse_iterator< ConstPartonIterator>; //! Default PhaseSpacePoint Constructor PhaseSpacePoint() = delete; //! 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. subl_chance gives * the chance of turning one emissions into a subleading configuration, * i.e. either unordered or central quark/anti-quark pair. Unordered * emissions require that the most extremal emission in any direction is * a quark or anti-quark and the next emission is a gluon. Quark/anti-quark * pairs are only generated for W processes. At most one subleading * emission will be generated in this way. */ PhaseSpacePoint( Process const & proc, JetParameters const & jet_properties, HEJ::PDF & pdf, double E_beam, double subl_chance, Subleading subl_channels, ParticlesDecayMap const & particle_decays, HEJ::EWConstants const & ew_parameters, 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_; } //! Iterator to the first outgoing parton ConstPartonIterator begin_partons() const; //! Iterator to the first outgoing parton ConstPartonIterator cbegin_partons() const; //! Iterator to the end of the outgoing partons ConstPartonIterator end_partons() const; //! Iterator to the end of the outgoing partons ConstPartonIterator cend_partons() const; //! Reverse Iterator to the first outgoing parton ConstReversePartonIterator rbegin_partons() const; //! Reverse Iterator to the first outgoing parton ConstReversePartonIterator crbegin_partons() const; //! Reverse Iterator to the first outgoing parton ConstReversePartonIterator rend_partons() const; //! Reverse Iterator to the first outgoing parton ConstReversePartonIterator crend_partons() const; private: friend HEJ::Event::EventData to_EventData(PhaseSpacePoint psp); /** * @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 ); HEJ::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( Process const & proc, double subl_chance, Subleading subl_channels, HEJ::PDF & pdf, double E_beam, double uf, HEJ::RNG & ran ); HEJ::ParticleID generate_incoming_id( std::size_t beam_idx, double x, double uf, HEJ::PDF & pdf, std::bitset<11> allowed_partons, HEJ::RNG & ran ); //! @brief Returns list of all allowed initial states partons std::array,2> filter_partons( Process const & proc, double const subl_chance, Subleading const subl_channels, HEJ::RNG & ran ); bool momentum_conserved(double ep) const; bool extremal_FKL_ok( std::vector const & partons ) const; double random_normal(double stddev, HEJ::RNG & ran); /** * @internal * @brief Turns a FKL configuration into a subleading one * * @param chance Chance to switch to subleading configuration * @param channels Allowed channels for subleading process * @param proc Process to decide which subleading * configurations are allowed * * With a chance of "chance" the FKL configuration is either turned into * a unordered configuration or, for A/W/Z bosons, a configuration with * a central quark/anti-quark pair. */ void maybe_turn_to_subl(double chance, Subleading channels, Process const & proc, HEJ::RNG & ran); void turn_to_uno(bool can_be_uno_backward, bool can_be_uno_forward, HEJ::RNG & ran); int select_qqx_flavour(bool allow_strange, HEJ::RNG & ran); - void turn_to_qqx(bool allow_strange, HEJ::RNG & ran); void turn_to_cqqx(bool allow_strange, HEJ::RNG & ran); void turn_to_eqqx(bool allow_strange, HEJ::RNG & ran); //! decay where we select the decay channel std::vector decay_boson( HEJ::Particle const & parent, std::vector const & decays, HEJ::RNG & ran ); //! generate decay products of a boson std::vector decay_boson( HEJ::Particle const & parent, std::vector const & decays, HEJ::RNG & ran ); /// @brief setup outgoing partons to ensure correct coupling to boson void couple_boson(HEJ::ParticleID boson, 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 ); //! Iterator over partons (non-const) using PartonIterator = boost::filter_iterator< bool (*)(HEJ::Particle const &), std::vector::iterator >; //! Reverse Iterator over partons (non-const) using ReversePartonIterator = std::reverse_iterator; //! Iterator to the first outgoing parton (non-const) PartonIterator begin_partons(); //! Iterator to the end of the outgoing partons (non-const) PartonIterator end_partons(); //! Reverse Iterator to the first outgoing parton (non-const) ReversePartonIterator rbegin_partons(); //! Reverse Iterator to the first outgoing parton (non-const) ReversePartonIterator rend_partons(); double weight_; Status status_; std::array incoming_; std::vector outgoing_; //! Particle decays in the format {outgoing index, decay products} std::unordered_map> decays_; }; //! Extract HEJ::Event::EventData from PhaseSpacePoint HEJ::Event::EventData to_EventData(PhaseSpacePoint psp); } diff --git a/FixedOrderGen/include/Subleading.hh b/FixedOrderGen/include/Subleading.hh index cb79ae5..8395370 100644 --- a/FixedOrderGen/include/Subleading.hh +++ b/FixedOrderGen/include/Subleading.hh @@ -1,30 +1,29 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019 * \copyright GPLv2 or later */ #pragma once #include namespace HEJFOG { namespace subleading { //!< @TODO confusing name with capital Subleading /** * Bit position of different subleading channels * e.g. (unsigned int) 1 => only unordered */ enum Channels: unsigned { - uno = 0u, + uno, unordered = uno, - qqx = 1u, - cqqx = 2u, + cqqx, central_qqx = cqqx, - eqqx = 3u, + eqqx, extremal_qqx = eqqx, first = uno, last = eqqx, }; } using Subleading = std::bitset; } diff --git a/FixedOrderGen/src/PhaseSpacePoint.cc b/FixedOrderGen/src/PhaseSpacePoint.cc index 4f75ab6..fad6bbe 100644 --- a/FixedOrderGen/src/PhaseSpacePoint.cc +++ b/FixedOrderGen/src/PhaseSpacePoint.cc @@ -1,940 +1,916 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "PhaseSpacePoint.hh" #include #include #include #include #include #include #include #include #include #include #include "fastjet/ClusterSequence.hh" #include "HEJ/Constants.hh" #include "HEJ/EWConstants.hh" #include "HEJ/exceptions.hh" #include "HEJ/kinematics.hh" #include "HEJ/Particle.hh" #include "HEJ/PDF.hh" #include "HEJ/RNG.hh" #include "HEJ/utility.hh" #include "JetParameters.hh" #include "Process.hh" namespace { using namespace HEJ; static_assert( std::numeric_limits::has_quiet_NaN, "no quiet NaN for double" ); constexpr double NaN = std::numeric_limits::quiet_NaN(); } // namespace anonymous namespace HEJFOG { Event::EventData to_EventData(PhaseSpacePoint psp){ //! @TODO Same function already in HEJ Event::EventData result; result.incoming = std::move(psp).incoming_; assert(result.incoming.size() == 2); result.outgoing = std::move(psp).outgoing_; // 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); result.decays = std::move(psp).decays_; result.parameters.central = {NaN, NaN, psp.weight()}; return result; } PhaseSpacePoint::ConstPartonIterator PhaseSpacePoint::begin_partons() const { return cbegin_partons(); } PhaseSpacePoint::ConstPartonIterator PhaseSpacePoint::cbegin_partons() const { return boost::make_filter_iterator( static_cast(HEJ::is_parton), cbegin(outgoing()), cend(outgoing()) ); } PhaseSpacePoint::ConstPartonIterator PhaseSpacePoint::end_partons() const { return cend_partons(); } PhaseSpacePoint::ConstPartonIterator PhaseSpacePoint::cend_partons() const { return boost::make_filter_iterator( static_cast(HEJ::is_parton), cend(outgoing()), cend(outgoing()) ); } PhaseSpacePoint::ConstReversePartonIterator PhaseSpacePoint::rbegin_partons() const { return crbegin_partons(); } PhaseSpacePoint::ConstReversePartonIterator PhaseSpacePoint::crbegin_partons() const { return std::reverse_iterator( cend_partons() ); } PhaseSpacePoint::ConstReversePartonIterator PhaseSpacePoint::rend_partons() const { return crend_partons(); } PhaseSpacePoint::ConstReversePartonIterator PhaseSpacePoint::crend_partons() const { return std::reverse_iterator( cbegin_partons() ); } PhaseSpacePoint::PartonIterator PhaseSpacePoint::begin_partons() { return boost::make_filter_iterator( static_cast(HEJ::is_parton), begin(outgoing_), end(outgoing_) ); } PhaseSpacePoint::PartonIterator PhaseSpacePoint::end_partons() { return boost::make_filter_iterator( static_cast(HEJ::is_parton), end(outgoing_), end(outgoing_) ); } PhaseSpacePoint::ReversePartonIterator PhaseSpacePoint::rbegin_partons() { return std::reverse_iterator( end_partons() ); } PhaseSpacePoint::ReversePartonIterator PhaseSpacePoint::rend_partons() { return std::reverse_iterator( begin_partons() ); } namespace { bool can_swap_to_uno( Particle const & p1, Particle const & p2 ){ assert(is_parton(p1) && is_parton(p2)); return p1.type != pid::gluon && p2.type == pid::gluon; } size_t count_gluons(PhaseSpacePoint::ConstPartonIterator first, PhaseSpacePoint::ConstPartonIterator last ){ return std::count_if(first, last, [](Particle const & p) {return p.type == pid::gluon;}); } /** assumes FKL configurations between first and last, * else there can be a quark in a non-extreme position * e.g. uno configuration gqg would pass */ Subleading possible_qqx( PhaseSpacePoint::ConstPartonIterator first, PhaseSpacePoint::ConstReversePartonIterator last ){ using namespace subleading; assert( std::distance( first,last.base() )>2 ); Subleading channels(~0l); - channels.reset(qqx); channels.reset(eqqx); channels.reset(cqqx); auto const ngluon = count_gluons(first,last.base()); if(ngluon < 2) return channels; - channels.set(qqx); if(first->type==pid::gluon || last->type==pid::gluon){ channels.set(eqqx); } if(std::distance(first,last.base())>=4){ channels.set(cqqx); } return channels; } bool is_AWZ_proccess(Process const & proc){ return proc.boson && is_AWZ_boson(*proc.boson); } bool is_up_type(Particle const & part){ return is_anyquark(part) && !(std::abs(part.type)%2); } bool is_down_type(Particle const & part){ return is_anyquark(part) && std::abs(part.type)%2; } bool can_couple_to_W(Particle const & part, pid::ParticleID const W_id){ const int W_charge = W_id>0?1:-1; return std::abs(part.type)<5 && ( (W_charge*part.type > 0 && is_up_type(part)) || (W_charge*part.type < 0 && is_down_type(part)) ); } } void PhaseSpacePoint::maybe_turn_to_subl( double chance, Subleading channels, Process const & proc, RNG & ran ){ if(proc.njets <= 2) return; assert(outgoing_.size() >= 2); // decide what kind of subleading process is allowed bool can_be_uno_backward = can_swap_to_uno( *cbegin_partons(), *(++cbegin_partons()) ); bool can_be_uno_forward = can_swap_to_uno( *crbegin_partons(), *(++crbegin_partons()) ); // Special case: Higgs can not be outside of uno if(proc.boson && *proc.boson==pid::Higgs){ if(outgoing_.begin()->type == pid::Higgs || (++outgoing_.begin())->type==pid::Higgs){ can_be_uno_backward = false; } if(outgoing_.rbegin()->type == pid::Higgs || (++outgoing_.rbegin())->type==pid::Higgs){ can_be_uno_forward = false; } } if(channels[subleading::uno]){ channels.set(subleading::uno, can_be_uno_backward || can_be_uno_forward); } channels &= possible_qqx(cbegin_partons(), crbegin_partons()); bool allow_strange = true; if(is_AWZ_proccess(proc)) { if(std::none_of(cbegin_partons(), cend_partons(), [&proc](Particle const & p){ return can_couple_to_W(p, *proc.boson);})) { // enforce qqx if A/W/Z can't couple somewhere else // this is ensured to work through filter_partons in reconstruct_incoming channels.reset(subleading::uno); assert(channels.any()); chance = 1.; // strange not allowed for W if(std::abs(*proc.boson)== pid::Wp) allow_strange = false; } } std::size_t const nchannels = channels.count(); // no subleading if(nchannels==0) return; if(ran.flat() >= chance){ weight_ /= 1 - chance; return; } weight_ /= chance; // select channel double const step = 1./nchannels; weight_*=nchannels; unsigned selected = subleading::first; double rnd = nchannels>1?ran.flat():0.; for(; selected<=subleading::last; ++selected){ assert(rnd>=0); if(channels[selected]){ if(rndtype, (++begin_partons())->type); } return std::swap(rbegin_partons()->type, (++rbegin_partons())->type); } if(can_be_uno_backward){ return std::swap(begin_partons()->type, (++begin_partons())->type); } assert(can_be_uno_forward); std::swap(rbegin_partons()->type, (++rbegin_partons())->type); } //! select flavour of quark int PhaseSpacePoint::select_qqx_flavour(const bool allow_strange, RNG & ran){ const double r1 = 2.*ran.flat()-1.; const double max_flavour = allow_strange?N_F:N_F-1; weight_ *= max_flavour*2; int flavour = pid::down + std::floor(std::abs(r1)*max_flavour); return flavour*(r1<0.?-1:1); } - void PhaseSpacePoint::turn_to_qqx(const bool allow_strange, RNG & ran){ - /// find first and last gluon in FKL chain - auto first = std::find_if(outgoing_.begin(), outgoing_.end(), - [](Particle const & p){return p.type == pid::gluon;}); - std::vector::iterator> FKL_gluons; - for(auto p = first; p!=outgoing_.end(); ++p){ - if(p->type == pid::gluon) FKL_gluons.push_back(p); - else if(is_anyquark(*p)) break; - } - const size_t ng = FKL_gluons.size(); - if(ng < 2) - throw std::logic_error("not enough gluons to create qqx"); - int flavour = select_qqx_flavour(allow_strange, ran); - // select gluon for switch - const size_t idx = std::floor((ng-1) * ran.flat()); - weight_ *= (ng-1); - FKL_gluons[idx]->type = ParticleID(flavour); - FKL_gluons[idx+1]->type = ParticleID(-flavour); - } - void PhaseSpacePoint::turn_to_cqqx(const bool allow_strange, RNG & ran){ // we assume all FKL partons to be gluons auto first = ++begin_partons(); auto last = ++rbegin_partons(); auto const ng = std::distance(first, last.base()); if(ng < 2) throw std::logic_error("not enough gluons to create qqx"); int flavour = select_qqx_flavour(allow_strange, ran); // select gluon for switch if(ng!=2){ const double steps = 1./(ng-1); weight_ /= steps; for(auto rnd = ran.flat(); rnd>steps; ++first){ rnd-=steps; } } first->type = static_cast(flavour); (++first)->type = static_cast(-flavour); } void PhaseSpacePoint::turn_to_eqqx(const bool allow_strange, RNG & ran){ /// find first and last gluon in FKL chain auto first = begin_partons(); const bool can_forward = !is_anyquark(*first); auto last = rbegin_partons(); const bool can_backward = !is_anyquark(*last); if(std::distance(first, last.base()) < 2) throw std::logic_error("not enough gluons to create qqx"); int flavour = select_qqx_flavour(allow_strange, ran); // select gluon for switch if(can_forward && !can_backward){ first->type = static_cast(flavour); (++first)->type = static_cast(-flavour); return; } if(!can_forward && can_backward){ last->type = static_cast(flavour); (++last)->type = static_cast(-flavour); return; } assert(can_forward && can_backward); weight_*=2.; if(ran.flat()>0.5){ first->type = static_cast(flavour); (++first)->type = static_cast(-flavour); return; } last->type = static_cast(flavour); (++last)->type = static_cast(-flavour); } 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 = std::sqrt(pt[0]*pt[0]+pt[1]*pt[1]+mass_square); const double pz=mperp*std::sinh(y); const double E=mperp*std::cosh(y); return {pt[0], pt[1], pz, E}; } namespace { //! adds a particle to target (in correct rapidity ordering) //! @returns positon of insertion auto insert_particle(std::vector & target, Particle && particle ){ const auto pos = std::upper_bound( begin(target),end(target),particle,rapidity_less{} ); target.insert(pos, std::move(particle)); return pos; } } PhaseSpacePoint::PhaseSpacePoint( Process const & proc, JetParameters const & jet_param, PDF & pdf, double E_beam, double const subl_chance, Subleading const subl_channels, ParticlesDecayMap const & particle_decays, EWConstants const & ew_parameters, RNG & ran ){ assert(proc.njets >= 2); status_ = Status::good; weight_ = 1; const int nout = proc.njets + (proc.boson?1:0) + proc.boson_decay.size(); outgoing_.reserve(nout); // generate parton momenta const bool is_pure_jets = (nout == proc.njets); auto partons = gen_LO_partons( proc.njets, is_pure_jets, jet_param, E_beam, ran ); // pre fill flavour with gluons for(auto&& p_out: partons) { outgoing_.emplace_back(Particle{pid::gluon, std::move(p_out), {}}); } if(status_ != Status::good) return; if(proc.boson){ // decay boson auto const & boson_prop = ew_parameters.prop(*proc.boson) ; auto boson{ gen_boson(*proc.boson, boson_prop.mass, boson_prop.width, ran) }; const auto pos{insert_particle(outgoing_, std::move(boson))}; const size_t boson_idx = std::distance(begin(outgoing_), pos); auto const & boson_decay = particle_decays.find(*proc.boson); if( !proc.boson_decay.empty() ){ // decay given in proc decays_.emplace( boson_idx, decay_boson(outgoing_[boson_idx], proc.boson_decay, ran) ); } else if( boson_decay != particle_decays.end() && !boson_decay->second.empty() ){ // decay given explicitly decays_.emplace( boson_idx, decay_boson(outgoing_[boson_idx], boson_decay->second, ran) ); } } // normalisation of momentum-conserving delta function weight_ *= std::pow(2*M_PI, 4); /** @TODO * uf (jet_param.min_pt) doesn'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 */ reconstruct_incoming(proc, subl_chance, subl_channels, pdf, E_beam, jet_param.min_pt, ran); if(status_ != Status::good) return; // set outgoing states begin_partons()->type = incoming_[0].type; rbegin_partons()->type = incoming_[1].type; maybe_turn_to_subl(subl_chance, subl_channels, proc, ran); if(proc.boson) couple_boson(*proc.boson, ran); } // pt generation, see eq:pt_sampling in developer manual double PhaseSpacePoint::gen_hard_pt( const int np , const double ptmin, const double ptmax, const double /* y */, RNG & ran ) { // heuristic parameter for pt sampling, see eq:pt_par in developer manual const double ptpar = ptmin + np/5.; const double arctan = std::atan((ptmax - ptmin)/ptpar); const double xpt = ran.flat(); const double pt = ptmin + ptpar*std::tan(xpt*arctan); const double cosine = std::cos(xpt*arctan); weight_ *= pt*ptpar*arctan/(cosine*cosine); return pt; } double PhaseSpacePoint::gen_soft_pt(int np, double max_pt, 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, 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_ = 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, RNG & ran ){ if (np<2) throw std::invalid_argument{"Not enough partons in gen_LO_partons"}; weight_ /= std::pow(16.*std::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_ = Status::not_enough_jets; return {}; } std::sort(begin(partons), end(partons), rapidity_less{}); return partons; } Particle PhaseSpacePoint::gen_boson( ParticleID bosonid, double mass, double width, RNG & ran ){ // Usual phase space measure weight_ /= 16.*std::pow(M_PI, 3); // Generate a y Gaussian distributed around 0 /// @TODO check magic numbers for different boson Higgs /// @TODO better sampling for W const double stddev_y = 1.6; const double y = random_normal(stddev_y, ran); const double r1 = ran.flat(); const double s_boson = mass*( mass + width*std::tan(M_PI/2.*r1 + (r1-1.)*std::atan(mass/width)) ); // off-shell s_boson sampling, compensates for Breit-Wigner /// @TODO use a flag instead if(std::abs(bosonid) == pid::Wp){ weight_/=M_PI*M_PI*16.; //Corrects B-W factors, see git issue 132 weight_*= mass*width*( M_PI+2.*std::atan(mass/width) ) / ( 1. + std::cos( M_PI*r1 + 2.*(r1-1.)*std::atan(mass/width) ) ); } auto p = gen_last_momentum(outgoing_, s_boson, y); return Particle{bosonid, std::move(p), {}}; } namespace { /// partons are ordered: even = anti, 0 = gluon ParticleID index_to_pid(size_t i){ if(!i) return pid::gluon; return static_cast( i%2 ? (i+1)/2 : -i/2 ); } /// partons are ordered: even = anti, 0 = gluon size_t pid_to_index(ParticleID id){ if(id==pid::gluon) return 0; return id>0 ? id*2-1 : std::abs(id)*2; } using part_mask = std::bitset<11>; //!< Selection mask for partons part_mask init_allowed(ParticleID const id){ if(std::abs(id) == pid::proton) return ~0; part_mask out{0}; if(is_parton(id)) out[pid_to_index(id)] = 1; return out; } /// decides which "index" (see index_to_pid) are allowed for process part_mask allowed_quarks(ParticleID const boson){ if(std::abs(boson) != pid::Wp){ return ~0; } // special case W: // Wp: anti-down or up-type quark, no b/t // Wm: down or anti-up-type quark, no b/t return boson>0?0b00011001101 :0b00100110011; } } /** * @brief Returns list of all allowed initial states partons * * checks which partons are allowed as initial state: * 1. only allow what is given in the Runcard (p -> all) * 2. A/W/Z require something to couple to * a) no qqx => no incoming gluon * b) 2j => no incoming gluon * c) 3j => can couple OR is gluon => 2 gluons become qqx later * 3. pure eqqx requires at least one gluon * 4. pure uno requires at least one quark */ std::array PhaseSpacePoint::filter_partons( Process const & proc, double const subl_chance, Subleading const subl_channels, RNG & ran ){ std::array allowed_partons{ init_allowed(proc.incoming[0]), init_allowed(proc.incoming[1]) }; // special case A/W/Z && central qqx not possible if(is_AWZ_proccess(proc) && ((proc.njets < 4) || !subl_channels[subleading::cqqx]) ){ // all possible incoming states auto allowed = allowed_quarks(*proc.boson); if(proc.njets == 2 || !subl_channels[subleading::eqqx]){ allowed.reset(0); } // possible states per leg std::array const maybe_partons{ allowed_partons[0]&allowed, allowed_partons[1]&allowed}; if(maybe_partons[0].any() && maybe_partons[1].any()){ // two options to get allowed initial state => choose one at random const size_t idx = ran.flat() < 0.5; allowed_partons[idx] = maybe_partons[idx]; //! @TODO Is this factor correct? I think so wasn't there before weight_*=2.; // else choose the possible } else if(maybe_partons[0].any()) { allowed_partons[0] = maybe_partons[0]; } else if(maybe_partons[1].any()) { allowed_partons[1] = maybe_partons[1]; } else{ throw std::invalid_argument{"Incoming state not allowed."}; } } // special case: pure subleading if(subl_chance!=1.){ return allowed_partons; } auto other_channels = subl_channels; // pure eqqx other_channels.reset(subleading::eqqx); if(other_channels.none()){ auto const gluon_idx = pid_to_index(pid::gluon); auto & first_beam = allowed_partons[0]; auto & second_beam = allowed_partons[1]; if(first_beam[gluon_idx] && !second_beam[gluon_idx]){ first_beam.reset(); first_beam.set(gluon_idx); return allowed_partons; } if(!first_beam[gluon_idx] && second_beam[gluon_idx]) { second_beam.reset(); second_beam.set(gluon_idx); return allowed_partons; } if(first_beam[gluon_idx] && second_beam[gluon_idx]) { // both beams can be gluons // if one can't be a quark everything is good auto first_quarks = first_beam; first_quarks.reset(gluon_idx); auto second_quarks = second_beam; second_quarks.reset(gluon_idx); if(first_quarks.none() || second_quarks.none()){ return allowed_partons; } // else choose one to be a gluon double rnd = ran.flat(); weight_*=3.; if(rnd<1./3.){ allowed_partons[0].reset(); allowed_partons[0].set(gluon_idx); allowed_partons[1].reset(gluon_idx); } else if(rnd<2./3.){ allowed_partons[1].reset(); allowed_partons[1].set(gluon_idx); allowed_partons[0].reset(gluon_idx); } else { allowed_partons[0].reset(); allowed_partons[0].set(gluon_idx); allowed_partons[1].reset(); allowed_partons[1].set(gluon_idx); } return allowed_partons; } throw std::invalid_argument{ "Incoming state not allowed for pure central qqx."}; } other_channels = subl_channels; // pure uno other_channels.reset(subleading::uno); if(other_channels.none()){ auto const gluon_idx = pid_to_index(pid::gluon); auto & first_beam = allowed_partons[0]; auto first_quarks = first_beam; first_quarks.reset(gluon_idx); auto & second_beam = allowed_partons[1]; auto second_quarks = second_beam; second_quarks.reset(gluon_idx); if(first_quarks.any() && second_quarks.none()){ first_beam.reset(gluon_idx); return allowed_partons; } if(first_quarks.none() && second_quarks.any()) { second_beam.reset(gluon_idx); return allowed_partons; } if(first_quarks.any() && second_quarks.any()) { // both beams can be quarks // if one can't be gluon everything is good if(!first_beam[gluon_idx] || !second_beam[gluon_idx]){ return allowed_partons; } // else choose one to be a quark double rnd = ran.flat(); weight_*=3.; if(rnd<1./3.){ allowed_partons[0].reset(gluon_idx); allowed_partons[1].reset(); allowed_partons[1].set(gluon_idx); } else if(rnd<2./3.){ allowed_partons[1].reset(gluon_idx); allowed_partons[0].reset(); allowed_partons[0].set(gluon_idx); } else { allowed_partons[0].reset(gluon_idx); allowed_partons[1].reset(gluon_idx); } return allowed_partons; } throw std::invalid_argument{ "Incoming state not allowed for pure unordered."}; } return allowed_partons; } void PhaseSpacePoint::reconstruct_incoming( Process const & proc, double const subl_chance, Subleading const subl_channels, PDF & pdf, double E_beam, double uf, 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].E()-incoming_[0].pz())/sqrts; const double xb=(incoming_[1].E()+incoming_[1].pz())/sqrts; // abort if phase space point is outside of collider energy reach if (xa>1. || xb>1.){ weight_=0; status_ = Status::too_much_energy; return; } auto const & ids = proc.incoming; std::array allowed_partons = filter_partons(proc, subl_chance, subl_channels, ran); for(size_t i = 0; i < 2; ++i){ if(ids[i] == pid::proton || ids[i] == pid::p_bar){ // pick ids according to pdfs incoming_[i].type = generate_incoming_id(i, i?xb:xa, uf, pdf, allowed_partons[i], ran); } else { assert(allowed_partons[i][pid_to_index(ids[i])]); incoming_[i].type = ids[i]; } } assert(momentum_conserved(1e-7)); } ParticleID PhaseSpacePoint::generate_incoming_id( size_t const beam_idx, double const x, double const uf, PDF & pdf, part_mask allowed_partons, RNG & ran ){ std::array pdf_wt; pdf_wt[0] = allowed_partons[0]? std::fabs(pdf.pdfpt(beam_idx,x,uf,pid::gluon)):0.; double pdftot = pdf_wt[0]; for(size_t i = 1; i < pdf_wt.size(); ++i){ pdf_wt[i] = allowed_partons[i]? 4./9.*std::fabs(pdf.pdfpt(beam_idx,x,uf,index_to_pid(i))):0; 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: "< allowed_parts; for(auto part_it=begin_partons(); part_it!=end_partons(); ++part_it){ // Wp -> up OR anti-down, Wm -> anti-up OR down, no bottom if ( can_couple_to_W(*part_it, boson) ) allowed_parts.push_back(part_it); } if(allowed_parts.size() == 0){ throw std::logic_error{"Found no parton for coupling with boson"}; } // select one and flip it size_t idx = 0; if(allowed_parts.size() > 1){ /// @TODO more efficient sampling /// old code: probability[i] = exp(parton[i].y - W.y) idx = std::floor(ran.flat()*allowed_parts.size()); weight_ *= allowed_parts.size(); } const int W_charge = boson>0?1:-1; allowed_parts[idx]->type = static_cast( allowed_parts[idx]->type - W_charge ); } double PhaseSpacePoint::random_normal( double stddev, RNG & ran ){ const double r1 = ran.flat(); const double r2 = ran.flat(); const double lninvr1 = -std::log(r1); const double result = stddev*std::sqrt(2.*lninvr1)*std::cos(2.*M_PI*r2); weight_ *= exp(result*result/(2*stddev*stddev))*std::sqrt(2.*M_PI)*stddev; return result; } 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); } Decay PhaseSpacePoint::select_decay_channel( std::vector const & decays, 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; if(decays.size()==1) return decays.front(); 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( Particle const & parent, std::vector const & decays, RNG & ran ){ const auto channel = select_decay_channel(decays, ran); return decay_boson(parent, channel.products, ran); } std::vector PhaseSpacePoint::decay_boson( Particle const & parent, std::vector const & decays, RNG & ran ){ if(decays.size() != 2){ throw not_implemented{ "only decays into two particles are implemented" }; } std::vector decay_products(decays.size()); for(size_t i = 0; i < decays.size(); ++i){ decay_products[i].type = decays[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.; // Jacobian Factors for W in line 418 const double sin_phi = std::sqrt(1. - cos_phi*cos_phi); // Know 0 < phi < pi const double px = E*std::cos(theta)*sin_phi; const double py = E*std::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/FixedOrderGen/src/config.cc b/FixedOrderGen/src/config.cc index 2d5be60..42c39f5 100644 --- a/FixedOrderGen/src/config.cc +++ b/FixedOrderGen/src/config.cc @@ -1,442 +1,440 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "config.hh" #include #include #include #include #include #include #include #include "HEJ/exceptions.hh" #include "HEJ/PDG_codes.hh" #include "HEJ/YAMLreader.hh" namespace HEJFOG { using HEJ::set_from_yaml; using HEJ::set_from_yaml_if_defined; using HEJ::pid::ParticleID; namespace { //! Get YAML tree of supported options /** * The configuration file is checked against this tree of options * in assert_all_options_known. */ YAML::Node const & get_supported_options(){ const static YAML::Node supported = [](){ YAML::Node supported; static const auto opts = { "process", "events", "subleading fraction","subleading channels", "scales", "scale factors", "max scale ratio", "pdf", "vev", "event output", "analyses", "analysis", "import scales" }; // add subnodes to "supported" - the assigned value is irrelevant for(auto && opt: opts) supported[opt] = ""; for(auto && jet_opt: {"min pt", "peak pt", "algorithm", "R", "max rapidity"}){ supported["jets"][jet_opt] = ""; } for(auto && particle_type: {"Higgs", "W", "Z"}){ for(auto && particle_opt: {"mass", "width"}){ supported["particle properties"][particle_type][particle_opt] = ""; } } for(auto && particle_type: {"Higgs", "Wp", "W+", "Wm", "W-", "Z"}){ for(auto && particle_opt: {"into", "branching ratio"}){ supported["decays"][particle_type][particle_opt] = ""; } } for(auto && opt: {"mt", "use impact factors", "include bottom", "mb"}){ supported["Higgs coupling"][opt] = ""; } for(auto && beam_opt: {"energy", "particles"}){ supported["beam"][beam_opt] = ""; } for(auto && unweight_opt: {"sample size", "max deviation"}){ supported["unweight"][unweight_opt] = ""; } for(auto && opt: {"name", "seed"}){ supported["random generator"][opt] = ""; } return supported; }(); return supported; } JetParameters get_jet_parameters( YAML::Node const & node, std::string const & entry ){ const auto p = HEJ::get_jet_parameters(node, entry); JetParameters result; result.def = p.def; result.min_pt = p.min_pt; set_from_yaml(result.max_y, node, entry, "max rapidity"); set_from_yaml_if_defined(result.peak_pt, node, entry, "peak pt"); if(result.peak_pt && *result.peak_pt <= result.min_pt) throw std::invalid_argument{ "Value of option 'peak pt' has to be larger than 'min pt'." }; return result; } Beam get_Beam( YAML::Node const & node, std::string const & entry ){ Beam beam; std::vector particles; set_from_yaml(beam.energy, node, entry, "energy"); set_from_yaml_if_defined(particles, node, entry, "particles"); if(! particles.empty()){ for(HEJ::ParticleID particle: particles){ if(particle != HEJ::pid::p && particle != HEJ::pid::p_bar){ throw std::invalid_argument{ "Unsupported value in option " + entry + ": particles:" " only proton ('p') and antiproton ('p_bar') beams are supported" }; } } if(particles.size() != 2){ throw std::invalid_argument{"Not exactly two beam particles"}; } beam.particles.front() = particles.front(); beam.particles.back() = particles.back(); } return beam; } std::vector split( std::string const & str, std::string const & delims ){ std::vector result; for(size_t begin, end = 0; end != str.npos;){ begin = str.find_first_not_of(delims, end); if(begin == str.npos) break; end = str.find_first_of(delims, begin + 1); result.emplace_back(str.substr(begin, end - begin)); } return result; } std::invalid_argument invalid_incoming(std::string const & what){ return std::invalid_argument{ "Incoming particle type " + what + " not supported," " incoming particles have to be 'p', 'p_bar' or partons" }; } std::invalid_argument invalid_outgoing(std::string const & what){ return std::invalid_argument{ "Outgoing particle type " + what + " not supported," " outgoing particles have to be 'j', 'photon', 'H', 'Wm', 'Wp', 'e-', 'e+', 'nu_e', 'nu_e_bar'" }; } HEJ::ParticleID reconstruct_boson_id( std::vector const & ids ){ assert(ids.size()==2); const int pidsum = ids[0] + ids[1]; if(pidsum == +1) { assert(HEJ::is_antilepton(ids[0])); if(HEJ::is_antineutrino(ids[0])) { throw HEJ::not_implemented{"lepton-flavour violating final state"}; } assert(HEJ::is_neutrino(ids[1])); // charged antilepton + neutrino means we had a W+ return HEJ::pid::Wp; } if(pidsum == -1) { assert(HEJ::is_antilepton(ids[0])); if(HEJ::is_neutrino(ids[1])) { throw HEJ::not_implemented{"lepton-flavour violating final state"}; } assert(HEJ::is_antineutrino(ids[0])); // charged lepton + antineutrino means we had a W- return HEJ::pid::Wm; } throw HEJ::not_implemented{ "final state with leptons "+name(ids[0])+" and "+name(ids[1]) +" not supported" }; } Process get_process( YAML::Node const & node, std::string const & entry ){ Process result; std::string process_string; set_from_yaml(process_string, node, entry); assert(! process_string.empty()); const auto particles = split(process_string, " \n\t\v=>"); if(particles.size() < 3){ throw std::invalid_argument{ "Bad format in option process: '" + process_string + "', expected format is 'in1 in2 => out1 ...'" }; } result.incoming.front() = HEJ::to_ParticleID(particles[0]); result.incoming.back() = HEJ::to_ParticleID(particles[1]); for(size_t i = 0; i < result.incoming.size(); ++i){ const HEJ::ParticleID in = result.incoming[i]; if( in != HEJ::pid::proton && in != HEJ::pid::p_bar && !HEJ::is_parton(in) ){ throw invalid_incoming(particles[i]); } } result.njets = 0; for(size_t i = result.incoming.size(); i < particles.size(); ++i){ assert(! particles[i].empty()); if(particles[i] == "j") ++result.njets; else if(std::isdigit(particles[i].front()) && particles[i].back() == 'j') result.njets += std::stoi(particles[i]); else{ const auto pid = HEJ::to_ParticleID(particles[i]); if(pid==HEJ::pid::Higgs || pid==HEJ::pid::Wp || pid==HEJ::pid::Wm){ if(result.boson) throw std::invalid_argument{ "More than one outgoing boson is not supported" }; if(!result.boson_decay.empty()) throw std::invalid_argument{ "Production of a boson together with a lepton is not supported" }; result.boson = pid; } else if (HEJ::is_anylepton(pid)){ // Do not accept more leptons, if two leptons are already mentioned if( result.boson_decay.size()>=2 ) throw std::invalid_argument{"Too many leptons provided"}; if(result.boson) throw std::invalid_argument{ "Production of a lepton together with a boson is not supported" }; result.boson_decay.emplace_back(pid); } else { throw invalid_outgoing(particles[i]); } } } if(result.njets < 2){ throw std::invalid_argument{ "Process has to include at least two jets ('j')" }; } if(!result.boson_decay.empty()){ std::sort(begin(result.boson_decay),end(result.boson_decay)); assert(!result.boson); result.boson = reconstruct_boson_id(result.boson_decay); } return result; } HEJFOG::Subleading to_subleading_channel(YAML::Node const & yaml){ std::string name; using namespace HEJFOG::subleading; Subleading channel; set_from_yaml(name, yaml); if(name == "none") return channel; if(name == "all") return channel.set(); if(name == "unordered" || name == "uno") return channel.set(uno); - if(name == "qqx") - return channel.set(qqx); - if(name == "cqqx") + if(name == "central qqx" || name == "cqqx") return channel.set(cqqx); - if(name == "eqqx") + if(name == "extremal qqx" || name == "eqqx") return channel.set(eqqx); throw HEJ::unknown_option("Unknown subleading channel '"+name+"'"); } Subleading get_subleading_channels(YAML::Node const & node){ using YAML::NodeType; using namespace HEJFOG::subleading; // all channels allowed by default if(!node) return ~0u; switch(node.Type()){ case NodeType::Undefined: return ~0u; case NodeType::Null: return 0u; case NodeType::Scalar: return to_subleading_channel(node); case NodeType::Map: throw HEJ::invalid_type{"map is not a valid option for subleading channels"}; case NodeType::Sequence: Subleading channels; for(auto && channel_node: node){ channels |= get_subleading_channels(channel_node); } return channels; } throw std::logic_error{"unreachable"}; } Decay get_decay(YAML::Node const & node, std::string const & entry, std::string const & boson ){ Decay decay; set_from_yaml(decay.products, node, entry, boson, "into"); decay.branching_ratio=1; set_from_yaml_if_defined(decay.branching_ratio, node, entry, boson, "branching ratio"); return decay; } std::vector get_decays(YAML::Node const & node, std::string const & entry, std::string const & boson ){ using YAML::NodeType; if(!node[entry][boson]) return {}; switch(node[entry][boson].Type()){ case NodeType::Null: case NodeType::Undefined: return {}; case NodeType::Scalar: throw HEJ::invalid_type{"value is not a list of decays"}; case NodeType::Map: return {get_decay(node, entry, boson)}; case NodeType::Sequence: std::vector result; for(auto && decay_str: node[entry][boson]){ result.emplace_back(get_decay(decay_str, entry, boson)); } return result; } throw std::logic_error{"unreachable"}; } ParticlesDecayMap get_all_decays(YAML::Node const & node, std::string const & entry ){ if(!node[entry]) return {}; if(!node[entry].IsMap()) throw HEJ::invalid_type{entry + " have to be a map"}; ParticlesDecayMap result; for(auto const & sub_node: node[entry]) { const auto boson = sub_node.first.as(); const auto id = HEJ::to_ParticleID(boson); result.emplace(id, get_decays(node, entry, boson)); } return result; } HEJ::ParticleProperties get_particle_properties( YAML::Node const & node, std::string const & entry, std::string const & boson ){ HEJ::ParticleProperties result; set_from_yaml(result.mass, node, entry, boson, "mass"); set_from_yaml(result.width, node, entry, boson, "width"); return result; } HEJ::EWConstants get_ew_parameters(YAML::Node const & node){ HEJ::EWConstants result; double vev; set_from_yaml(vev, node, "vev"); result.set_vevWZH(vev, get_particle_properties(node, "particle properties", "W"), get_particle_properties(node, "particle properties", "Z"), get_particle_properties(node, "particle properties", "Higgs") ); return result; } UnweightSettings get_unweight( YAML::Node const & node, std::string const & entry ){ UnweightSettings result; set_from_yaml(result.sample_size, node, entry, "sample size"); if(result.sample_size <= 0){ throw std::invalid_argument{ "negative sample size " + std::to_string(result.sample_size) }; } set_from_yaml(result.max_dev, node, entry, "max deviation"); return result; } Config to_Config(YAML::Node const & yaml){ try{ HEJ::assert_all_options_known(yaml, get_supported_options()); } catch(HEJ::unknown_option const & ex){ throw HEJ::unknown_option{std::string{"Unknown option '"} + ex.what() + "'"}; } Config config; config.process = get_process(yaml, "process"); set_from_yaml(config.events, yaml, "events"); config.jets = get_jet_parameters(yaml, "jets"); config.beam = get_Beam(yaml, "beam"); for(size_t i = 0; i < config.process.incoming.size(); ++i){ auto const & in = config.process.incoming[i]; using namespace HEJ::pid; if( (in == p || in == p_bar) && in != config.beam.particles[i]){ throw std::invalid_argument{ "Particle type of beam " + std::to_string(i+1) + " incompatible" + " with type of incoming particle " + std::to_string(i+1) }; } } set_from_yaml(config.pdf_id, yaml, "pdf"); set_from_yaml(config.subleading_fraction, yaml, "subleading fraction"); if(config.subleading_fraction == 0) config.subleading_channels.reset(); else config.subleading_channels = get_subleading_channels(yaml["subleading channels"]); config.ew_parameters = get_ew_parameters(yaml); config.particle_decays = get_all_decays(yaml, "decays"); if(config.process.boson // check that Ws always decay && std::abs(*config.process.boson) == HEJ::ParticleID::Wp && config.process.boson_decay.empty() ){ auto const & decay = config.particle_decays.find(*config.process.boson); if(decay == config.particle_decays.end() || decay->second.empty()) throw std::invalid_argument{ "Decay for "+name(*config.process.boson)+" is required"}; } set_from_yaml_if_defined(config.analyses_parameters, yaml, "analyses"); if(yaml["analysis"]){ std::cerr << "WARNING: Configuration entry 'analysis' is deprecated. " " Use 'analyses' instead.\n"; YAML::Node ana; set_from_yaml(ana, yaml, "analysis"); if(!ana.IsNull()){ config.analyses_parameters.push_back(ana); } } config.scales = HEJ::to_ScaleConfig(yaml); set_from_yaml_if_defined(config.output, yaml, "event output"); config.rng = HEJ::to_RNGConfig(yaml, "random generator"); config.Higgs_coupling = HEJ::get_Higgs_coupling(yaml, "Higgs coupling"); if(yaml["unweight"]) config.unweight = get_unweight(yaml, "unweight"); return config; } } // namespace anonymous Config load_config(std::string const & config_file){ try{ return to_Config(YAML::LoadFile(config_file)); } catch(...){ std::cerr << "Error reading " << config_file << ":\n "; throw; } } } diff --git a/FixedOrderGen/t/W_nj_classify.cc b/FixedOrderGen/t/W_nj_classify.cc index a27225c..b2e02dd 100644 --- a/FixedOrderGen/t/W_nj_classify.cc +++ b/FixedOrderGen/t/W_nj_classify.cc @@ -1,228 +1,227 @@ /** * \brief check that the PSP generates the all W+jet subleading processes * * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include #include #include #include #include #include #include #include #include "HEJ/Event.hh" #include "HEJ/event_types.hh" #include "HEJ/EWConstants.hh" #include "HEJ/exceptions.hh" #include "HEJ/Mixmax.hh" #include "HEJ/PDF.hh" #include "HEJ/PDG_codes.hh" #include "fastjet/JetDefinition.hh" #include "Decay.hh" #include "JetParameters.hh" #include "PhaseSpacePoint.hh" #include "Process.hh" #include "Status.hh" #include "Subleading.hh" namespace { using namespace HEJFOG; using namespace HEJ; void print_psp(PhaseSpacePoint const & psp){ std::cerr << "Process:\n" << psp.incoming()[0].type << " + "<< psp.incoming()[1].type << " -> "; for(auto const & out: psp.outgoing()){ std::cerr << out.type << " "; } std::cerr << "\n"; } void bail_out(PhaseSpacePoint const & psp, std::string msg){ print_psp(psp); throw std::logic_error{msg}; } #if !defined(__clang__) && defined(__GNUC__) && (__GNUC__ < 6) // gcc version < 6 explicitly needs hash function for enum // see https://gcc.gnu.org/bugzilla/show_bug.cgi?id=60970 using type_map = std::unordered_map>; #else using type_map = std::unordered_map; #endif } int main(){ constexpr std::size_t n_psp_base = 10375; const JetParameters jet_para{ fastjet::JetDefinition(fastjet::JetAlgorithm::antikt_algorithm, 0.4), 30, 5, 30}; PDF pdf(11000, pid::proton, pid::proton); constexpr double E_cms = 13000.; constexpr double subl_chance = 0.8; const ParticlesDecayMap boson_decays{ {pid::Wp, {Decay{ {pid::e_bar, pid::nu_e}, 1.} }}, {pid::Wm, {Decay{ {pid::e, pid::nu_e_bar}, 1.} }} }; const EWConstants ew_constants{246.2196508, ParticleProperties{80.385, 2.085}, ParticleProperties{91.187, 2.495}, ParticleProperties{125, 0.004165} }; HEJ::Mixmax ran{}; Subleading subl_channels(~0l); std::vector allowed_types{event_type::FKL, event_type::unob, event_type::unof, event_type::qqxexb, event_type::qqxexf}; std::cout << "Wp3j" << std::endl; // Wp3j Process proc {{pid::proton,pid::proton}, 3, pid::Wp, {}}; std::size_t n_psp = n_psp_base; type_map type_counter; for( std::size_t i = 0; i try again ++n_psp; } } std::cout << "Wp+3j: Took " << n_psp << " to generate " << n_psp_base << " successfully PSP (" << 1.*n_psp/n_psp_base << " trials/PSP)" << std::endl; std::cout << "States by classification:\n"; for(auto const & entry: type_counter){ const double fraction = static_cast(entry.second)/n_psp_base; const int percent = std::round(100*fraction); std::cout << std::left << std::setw(25) << (name(entry.first) + std::string(":")) << entry.second << " (" << percent << "%)\n"; } for(auto const & t: allowed_types){ if(type_counter[t] < 0.05 * n_psp_base){ std::cerr << "Less than 5% of the events are of type " << name(t) << std::endl; return EXIT_FAILURE; } } // Wm3j - only uno proc = Process{{pid::proton,pid::proton}, 3, pid::Wm, {}}; n_psp = n_psp_base; subl_channels.reset(); subl_channels.set(subleading::uno); allowed_types = {event_type::FKL, event_type::unob, event_type::unof}; type_counter.clear(); for( std::size_t i = 0; i try again ++n_psp; } } std::cout << "Wm+3j (only uno): Took " << n_psp << " to generate " << n_psp_base << " successfully PSP (" << 1.*n_psp/n_psp_base << " trials/PSP)" << std::endl; std::cout << "States by classification:\n"; for(auto const & entry: type_counter){ const double fraction = static_cast(entry.second)/n_psp_base; const int percent = std::round(100*fraction); std::cout << std::left << std::setw(25) << (name(entry.first) + std::string(":")) << entry.second << " (" << percent << "%)\n"; } for(auto const & t: allowed_types){ if(type_counter[t] < 0.05 * n_psp_base){ std::cerr << "Less than 5% of the events are of type " << name(t) << std::endl; return EXIT_FAILURE; } } // Wm4j proc = Process{{pid::proton,pid::proton}, 4, pid::Wm, {}}; n_psp = n_psp_base; subl_channels.set(); - subl_channels.reset(subleading::qqx); //!< @TODO remove, temp fix to not double count allowed_types = {event_type::FKL, event_type::unob, event_type::unof, event_type::qqxexb, event_type::qqxexf, event_type::qqxmid}; type_counter.clear(); std::array wpos_counter; // position of the W boson (back, central, forward) for( std::size_t i = 0; itype==pid::Wm){ ++wpos_counter[0][ev.type()]; } else if(ev.outgoing().rbegin()->type==pid::Wm){ ++wpos_counter[2][ev.type()]; } else { ++wpos_counter[1][ev.type()]; } } else { // bad process -> try again ++n_psp; } } std::cout << "Wm+4j: Took " << n_psp << " to generate " << n_psp_base << " successfully PSP (" << 1.*n_psp/n_psp_base << " trials/PSP)" << std::endl; std::cout << "States by classification:\n"; for(auto const & entry: type_counter){ const double fraction = static_cast(entry.second)/n_psp_base; const int percent = std::round(100*fraction); std::cout << std::left << std::setw(25) << (name(entry.first) + std::string(":")) << entry.second << " (" << percent << "%)\n"; } for(auto const & t: allowed_types){ if(type_counter[t] < 0.03 * n_psp_base){ std::cerr << "Less than 3% of the events are of type " << name(t) << std::endl; return EXIT_FAILURE; } } std::cout << "Stats by Wpos:\n"; for(std::size_t i=0; i 2j subleading fraction: 0.1 subleading channels: - - eqqx - - uno + - extremal qqx + - central qqx + - unordered scales: 125 random generator: name: mixmax particle properties: Higgs: mass: 125 width: 0 W: mass: 80.385 width: 2.085 Z: mass: 91.187 width: 2.495 vev: 246.2196508 unweight: sample size: 100 max deviation: -10 diff --git a/FixedOrderGen/t/config_2j_peak.yml b/FixedOrderGen/t/config_2j_peak.yml index 346e13f..deeec6e 100644 --- a/FixedOrderGen/t/config_2j_peak.yml +++ b/FixedOrderGen/t/config_2j_peak.yml @@ -1,41 +1,41 @@ events: 1000 jets: min pt: 20 R: 0.4 algorithm: antikt max rapidity: 5 peak pt: 30 unweight: sample size: 100 max deviation: 1 beam: energy: 6500 particles: [p, p] pdf: 11000 process: p p => 2j subleading fraction: 0.1 -subleading channels: none +subleading channels: all scales: 125 random generator: name: mixmax particle properties: Higgs: mass: 125 width: 0 W: mass: 80.385 width: 2.085 Z: mass: 91.187 width: 2.495 vev: 246.2196508 diff --git a/doc/sphinx/HEJFOG.rst b/doc/sphinx/HEJFOG.rst index 7352786..b485fe4 100644 --- a/doc/sphinx/HEJFOG.rst +++ b/doc/sphinx/HEJFOG.rst @@ -1,296 +1,299 @@ The HEJ Fixed Order Generator ============================= For high jet multiplicities event generation with standard fixed-order generators becomes increasingly cumbersome. For example, the leading-order production of a Higgs Boson with five or more jets is computationally prohibitively expensive. To this end, HEJ 2 provides the ``HEJFOG`` fixed-order generator that allows to generate events with high jet multiplicities. To facilitate the computation the limit of Multi-Regge Kinematics with large invariant masses between all outgoing particles is assumed in the matrix elements. The typical use of the ``HEJFOG`` is to supplement low-multiplicity events from standard generators with high-multiplicity events before using the HEJ 2 program to add high-energy resummation. Installation ------------ The ``HEJFOG`` comes bundled together with HEJ 2 and the installation is very similar. After downloading HEJ 2 and installing the prerequisites as described in :ref:`Installation` the ``HEJFOG`` can be installed with:: cmake /path/to/FixedOrderGen -DCMAKE_INSTALL_PREFIX=target/directory make install where :file:`/path/to/FixedOrderGen` refers to the :file:`FixedOrderGen` subdirectory in the HEJ 2 directory. If HEJ 2 was installed to a non-standard location, it may be necessary to specify the directory containing :file:`HEJ-config.cmake`. If the base installation directory is :file:`/path/to/HEJ`, :file:`HEJ-config.cmake` should be found in :file:`/path/to/HEJ/lib/cmake/HEJ` and the commands for installing the ``HEJFOG`` would read:: cmake /path/to/FixedOrderGen -DHEJ_DIR=/path/to/HEJ/lib/cmake/HEJ -DCMAKE_INSTALL_PREFIX=target/directory make install The installation can be tested with:: make test provided that the ``CT10nlo`` PDF sets is installed. Running the fixed-order generator --------------------------------- After installing the ``HEJFOG`` you can modify the provided configuration file :file:`configFO.yml` and run the generator with:: HEJFOG configFO.yml The resulting event file, by default named :file:`HEJFO.lhe`, can then be fed into HEJ 2 like any event file generated from a standard fixed-order generator, see :ref:`Running HEJ 2`. Settings -------- Similar to HEJ 2, the ``HEJFOG`` uses a `YAML `_ configuration file. The settings are .. _`process`: **process** The scattering process for which events are being generated. The format is :code:`in1 in2 => out1 out2 ...` The incoming particles, :code:`in1`, :code:`in2` can be - quarks: :code:`u`, :code:`d`, :code:`u_bar`, and so on - gluons: :code:`g` - protons :code:`p` or antiprotons :code:`p_bar` The outgoing particles, can be - jets: :code:`j`, multiple jets can be grouped together e.g. :code:`2j` - bosons: Any gauge boson, they might further :ref:`decay` - leptons: if they can be reconstructed to a :code:`W+` or :code:`W-`, e.g. :code:`e+ nu_e` At most one of the outgoing particles can be a boson, the rest has to be partonic. At the moment only Higgs :code:`h`, :code:`W+` and :code:`W-` bosons are supported. All other outgoing particles are jets. There have to be at least two jets. Instead of the leptons decays of the bosons can be added through the :ref:`decays`. The latter is required to decay a Higgs boson. So :code:`p p => Wp j j` is the same as :code:`p p => e+ nu_e 2j`, if the decay :code:`Wp => e+ nu_e` is specified in `decays`_. .. _`events`: **events** Specifies the number of events to generate. .. _`jets`: **jets** Defines the properties of the generated jets. .. _`jets: min pt`: **min pt** Minimum jet transverse momentum in GeV. .. _`jets: peak pt`: **peak pt** Optional setting to specify the dominant jet transverse momentum in GeV. If the generated events are used as input for HEJ resummation, this should be set to the minimum transverse momentum of the resummation jets. In all cases it has to be larger than :code:`min pt`. The effect is that only a small fraction of jets will be generated with a transverse momentum below the value of this setting. .. _`jets: algorithm`: **algorithm** The algorithm used to define jets. Allowed settings are :code:`kt`, :code:`cambridge`, :code:`antikt`, :code:`cambridge for passive`. See the `FastJet `_ documentation for a description of these algorithms. .. _`jets: R`: **R** The R parameter used in the jet algorithm. .. _`jets: max rapidity`: **max rapidity** Maximum absolute value of the jet rapidity. .. _`beam`: **beam** Defines various properties of the collider beam. .. _`beam: energy`: **energy** The beam energy in GeV. For example, the 13 TeV LHC corresponds to a value of 6500. .. _`beam: particles`: **particles** A list :code:`[p1, p2]` of two beam particles. The only supported entries are protons :code:`p` and antiprotons :code:`p_bar`. .. _`pdf`: **pdf** The `LHAPDF number `_ of the PDF set. For example, 230000 corresponds to an NNPDF 2.3 NLO PDF set. .. _`subleading fraction`: **subleading fraction** This setting is related to the fraction of events that are not a FKL configuration. All possible subleading process are listed in :ref:`subleading channels`. This value must be positive and not larger than 1. It should typically be chosen between 0.01 and 0.1. Note that while this parameter influences the chance of generating subleading configurations, it generally does not correspond to the actual fraction of subleading events. .. _`subleading channels`: **subleading channels** Optional parameters to select a specific subleading process, multiple channels can be selected at once. If multiple subleading configurations are possible one will be selected at random for each event, thus each final state will include at most one subleading process, e.g. either - :code:`unordered` or :code:`qqx`. Only has an effect if + :code:`unordered` or :code:`central qqx`. Only has an effect if :code:`subleading fraction` is non-zero. Subleading processes requested here also have to be available from HEJ, for a list of implemented processes see :ref:`event treatment`. The following values are allowed: - :code:`all`: All channels allowed. This is the default. - :code:`none`: No subleading contribution, only the leading (FKL) configurations are allowed. This is equivalent to :code:`subleading fraction: 0`. - :code:`unordered`: Unordered emission allowed. Unordered emission are any rapidity ordering where exactly one gluon is emitted outside the rapidity ordering required in FKL events. More precisely, if at least one of the incoming particles is a quark or antiquark and there are more than two jets in the final state, :code:`subleading fraction` states the probability that the flavours of the outgoing particles are assigned in such a way that an unordered configuration arises. - - :code:`qqx`: Production of a quark-antiquark pair. In the leading - (FKL) configurations all partons except for the most backward and - forward ones are gluons. If the :code:`qqx` channel is allowed, + - :code:`extremal qqx`: Production of a quark-antiquark pair as extremal + partons in rapidity. If the :code:`central qqx` channel is allowed, :code:`subleading fraction` gives the probability of emitting a - quark-antiquark pair in place of two gluons that are adjacent in - rapidity. + quark-antiquark pair in place of two gluons adjacent in rapidity at the most forward or backwards quarks. + - :code:`central qqx`: Production of a central quark-antiquark pair. This + setting is similar to :code:`extremal qqx`, but :code:`subleading fraction` + gives the probability of emitting a quark-antiquark pair in place of two + gluons adjacent in rapidity at any point *inside* the :code:`FKL` gluon + chain. .. _`unweight`: **unweight** This setting defines the parameters for the partial unweighting of events. You can disable unweighting by removing this entry from the configuration file. .. _`unweight: sample size`: **sample size** The number of weighted events used to calibrate the unweighting. A good default is to set this to the number of target `events`_. If the number of `events`_ is large this can lead to significant memory consumption and a lower value should be chosen. Contrarily, for large multiplicities the unweighting efficiency becomes worse and the sample size should be increased. .. _`unweight: max deviation`: **max deviation** Controls the range of events to which unweighting is applied. A larger value means that a larger fraction of events are unweighted. Typical values are between -1 and 1. .. _`decays`: **decays** Optional setting specifying the decays of the particle. Only the decay into two particles is implemented. Each decay has the form :code:`boson: {into: [p1,p2], branching ratio: r}` where :code:`boson` is the corresponding boson, :code:`p1` and :code:`p2` are the particle names of the decay product (e.g. :code:`photon`) and :code:`r` is the branching ratio. The branching ratio is optional (:code:`1` by default). Decays of a Higgs boson are treated as the production and subsequent decay of an on-shell Higgs boson, so decays into e.g. Z bosons are not supported. In contrast W bosons are decayed off-shell, thus the branching ratio should not be specified or set to :code:`1`. .. _`FOG scales`: **scales** Specifies the renormalisation and factorisation scales for the output events. For details, see the corresponding entry in the HEJ 2 :ref:`settings`. Note that this should usually be a single value, as the weights resulting from additional scale choices will simply be ignored when adding high-energy resummation with HEJ 2. .. _`FOG event output`: **event output** Specifies the name of a single event output file or a list of such files. See the corresponding entry in the HEJ 2 :ref:`settings` for details. .. _`RanLux init`: .. _`FOG random generator`: **random generator** Sets parameters for random number generation. See the HEJ 2 :ref:`settings` for details. .. _`FOG analyses`: **analyses** Specifies the name and settings for one or more analyses library. This can be useful to specify cuts at the fixed-order level. See the corresponding entry in the HEJ 2 :ref:`settings` for details. .. _`FOG vev`: **vev** Higgs vacuum expectation value in GeV. All electro-weak constants are derived from this together with the :ref:`particle properties `. .. _`FOG particle properties`: **particle properties** Specifies various properties of the different particles (Higgs, W or Z). These have to be specified for all bosons. See the corresponding entry in the HEJ 2 :ref:`settings` for details. .. _`FOG Higgs coupling`: **Higgs coupling** This collects a number of settings concerning the effective coupling of the Higgs boson to gluons. See the corresponding entry in the HEJ 2 :ref:`settings` for details.