Page MenuHomeHEPForge
Files F8725507

No OneTemporary
Actions

View Options

diff --git a/include/RHEJ/CombinedEventWriter.hh b/include/RHEJ/CombinedEventWriter.hh
index 5f7c471..44308a0 100644
--- a/include/RHEJ/CombinedEventWriter.hh
+++ b/include/RHEJ/CombinedEventWriter.hh
@@ -1,53 +1,53 @@
/** \file CombinedEventWriter.hh
* \brief Details the CombinedEventWriter class
*
* CombinedEventWriter Class which handles all of the EventWriters being used by a single use
* of RHEJ. It calls all other EventWriters in its list to use their write function and output.
*/
#pragma once
#include <memory>
#include <vector>
-#include "RHEJ/config.hh"
#include "RHEJ/EventWriter.hh"
#include "RHEJ/make_writer.hh"
namespace LHEF{
/** \struct HEPRUP CombinedEventWriter.hh "include/RHEJ/CombinedEventWriter.hh"
* \brief HEPRUP struct which contains input event information.
*
* This struct contains important information from the LHEF header such as the beam energy.
*/
struct HEPRUP;
}
namespace RHEJ{
+ struct config;
/** \class CombinedEventWriter CombinedEventWriter.hh "include/RHEJ/CombinedEventWriter.hh"
* \brief Intended as an EventWriter to Multiple Output types.
*
* Inherits from EventWriter, and then uses its write function which calls all of the event
* writers write function to output to their specific output files in their specific format.
* This handles this complexity of having multiple EventWriters.
*/
class CombinedEventWriter: public EventWriter{
public:
//! CombinedEventWriter Constructor
CombinedEventWriter(
/** Contains all configurations for running rHEJ */
Config const & conf,
/** Contains important information from the Input */
LHEF::HEPRUP const & heprup
);
void write(Event const &) override;
private:
//! A vector of the EventWriters
std::vector<std::unique_ptr<EventWriter>> writers_;
};
}
diff --git a/include/RHEJ/EventReweighter.hh b/include/RHEJ/EventReweighter.hh
index 1463f21..22229f8 100644
--- a/include/RHEJ/EventReweighter.hh
+++ b/include/RHEJ/EventReweighter.hh
@@ -1,181 +1,175 @@
/** \file EventReweighter.hh
* \brief Main class for reweighting events according to the reversed HEJ method
*
* EventReweighter Class is the main class used within RHEJ. It reweights the
* resummation events.
*/
#pragma once
#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequence.hh"
#include "LHEF/LHEF.h"
+#include "RHEJ/config.hh"
#include "RHEJ/PDF.hh"
#include "RHEJ/utility.hh"
#include "RHEJ/Event.hh"
#include "RHEJ/event_types.hh"
-#include "RHEJ/config.hh"
#include "RHEJ/PhaseSpacePoint.hh"
#include "RHEJ/MatrixElement.hh"
namespace RHEJ{
- class ScaleGenerator;
/** \struct Beam EventReweighter.hh "include/RHEJ/EventReweighter.hh"
* \brief Beam Struct. Contains beam energy and incoming particle content.
*
* Contains the value of the energy of the beam and the identity of the two incoming Particles
*/
struct Beam{
double E;
std::array<ParticleID, 2> type;
};
/** \class EventReweighter EventReweighter.hh "include/RHEJ/EventReweighter.hh"
* \brief EventReweighter Main class for reweighting events in RHEJ.
*
* This is the class which reweights the resummation phase space points once they have been
* generated by PhaseSpacePoint. This class also handles how many phase space points are generated
* for every fixed order event which is processed.
*/
class EventReweighter{
using EventType = event_type::EventType;
public:
EventReweighter(
Beam beam, /**< Beam Energy */
int pdf_id, /**< PDF ID */
Config const & conf /**< Configuration parameters */
);
EventReweighter(
LHEF::HEPRUP const & heprup, /**< LHEF event header */
Config const & conf /**< Configuration parameters */
);
PDF const & pdf() const; /**< PDF Used */
std::vector<Event> reweight(
Event const & ev, /**< Event For Reweighing */
int num_events /**< Number of Events */
);
std::vector<Event> reweight(
Event const & ev, /**< Event For Reweighing */
int num_events, /**< Number of Events */
ScaleGenerator const & gen_scales /**< Generated Scales see config.hh */
);
private:
/** \struct EventFactors EventReweighter.hh "include/RHEJ/EventReweighter.hh"
* \brief EventFactors Struct containing central value and variations
*
* Contains the Central Value and Variation around
*/
struct EventFactors{
double central; /**< Central Value for the Factor */
std::vector<double> variations; /**< Vector of Variations from central value */
};
template<typename... T>
PDF const & pdf(T&& ...);
std::vector<Event> gen_res_events(
Event const & ev, int num_events
);
std::vector<Event> rescale(
Event const & Born_ev, std::vector<Event> events
) const;
/**
* \brief Do the Jets pass the resummation Cuts?
*
* @param ev Event in Question
* @returns 0 or 1 depending on if ev passes Jet Cuts
*/
bool jets_pass_resummation_cuts(Event const & ev) const;
/**
* \brief pdf_factors Function
*
* @param ev Event in Question
* @returns EventFactor due to PDFs
*
* Calculates the Central value and the variation due
* to the PDF choice made.
*/
EventFactors pdf_factors(Event const & ev) const;
/**
* \brief matrix_elements Function
*
* @param ev Event in question
* @returns EventFactor due to MatrixElements
*
* Calculates the Central value and the variation due
* to the Matrix Element.
*/
EventFactors matrix_elements(Event const & ev) const;
/**
* \brief Scale-dependent part of fixed-order matrix element
*
* @param ev Event in question
* @returns EventFactor scale variation due to FO-ME.
*
* This is only called to compute the scale variation for events where
* we don't do resummation (e.g. non-FKL).
* Since at tree level the scale dependence is just due to alpha_s,
* it is enough to return the alpha_s(mur) factors in the matrix element.
* The rest drops out in the ratio of (output event ME)/(input event ME),
* so we never have to compute it.
*/
EventFactors fixed_order_scale_ME(Event const & ev) const;
/**
* \brief Computes the tree level matrix element
*
* @param ev Event in Question
* @returns HEJ approximation to Tree level Matrix Element
*
* This computes the HEJ approximation to the tree level FO
* Matrix element which is used within the LO weighting process.
*/
double tree_matrix_element(Event const & ev) const;
/**
* \brief Generation of Resummation Phase Space Points
*
* @param ev Event In Question
* @returns Resummation PhaseSpacePoint
*
* Where is this even called? I can't find it with a grep -r in top directory???
*/
PhaseSpacePoint gen_PhaseSpacePoint(Event const & ev) const;
- const EventReweighterParameters param;
+ EventReweighterParameters param_;
- // const double extpartonptmin_; /**< External Parton Minimum Pt */
- // const double max_ext_soft_pt_fraction_; /**< External Parton Maximum Pt fraction */
- // const double jetptmin_; /**< Jet Pt Minimum threshold */
- const double E_beam_; /**< Energy of Beam */
+ double E_beam_; /**< Energy of Beam */
- // const fastjet::JetDefinition jet_def_; /**< Jet Definition being used */
- const PDF pdf_; /**< Relevant PDF */
+ PDF pdf_; /**< Relevant PDF */
- const MatrixElement MEt2_; /**< Matrix Element Object */
- // const EventTreatMap treat_;
+ MatrixElement MEt2_; /**< Matrix Element Object */
};
template<typename... T>
PDF const & EventReweighter::pdf(T&&... t){
return pdf_ = PDF{std::forward<T>(t)...};
}
}
diff --git a/include/RHEJ/EventWriter.hh b/include/RHEJ/EventWriter.hh
index e8f969c..c121049 100644
--- a/include/RHEJ/EventWriter.hh
+++ b/include/RHEJ/EventWriter.hh
@@ -1,29 +1,24 @@
/** \file EventWriter.hh
* \brief Contains the EventWriter Class which is inherited by all writers.
*/
#pragma once
#include <string>
#include <exception>
namespace RHEJ{
class Event;
- struct Config;
- /** \class EventWriter EventWriter.hh "include/RHEJ/EventWriter.hh"
+ /** \struct EventWriter EventWriter.hh "include/RHEJ/EventWriter.hh"
* \brief The basis of every other writer in RHEJ
*
* All other writers inherit from this, regardless of final
* output. Includes the virtual function write and a destructor
*/
- class EventWriter{
- protected:
- const Config & config;
- public:
- EventWriter(Config const & conf): config(conf){}
+ struct EventWriter{
virtual void write(Event const &) = 0;
virtual ~EventWriter() = default;
};
}
diff --git a/include/RHEJ/HepMCWriter.hh b/include/RHEJ/HepMCWriter.hh
index 8a46b7f..c6d2f0d 100644
--- a/include/RHEJ/HepMCWriter.hh
+++ b/include/RHEJ/HepMCWriter.hh
@@ -1,47 +1,47 @@
/** \file HepMCWriter.hh
* \brief Contains the EventWriter which is necessary for HepMC Output.
*/
#pragma once
#include <string>
-#include "RHEJ/config.hh"
#include "RHEJ/EventWriter.hh"
#include "LHEF/LHEF.h"
namespace RHEJ{
class Event;
+ struct Config;
/** \class HepMCWriter HepMCWriter.hh "include/RHEJ/HepMCWriter.hh"
* \brief This is an event writer specifically for HepMC output.
* as such it inherits everything from the EventWriter class
* and also includes HepMCWriter constructors
*
* Implementation note:
* This uses the pimpl ("pointer to implementation") idiom.
* HepMC support is optional and the implementation depends on the
* HepMC version. Without pimpl, we would have to specify the HepMC version
* via the preprocessor whenever this header is included. We don't want to
* burden users of the rHEJ library (for example the HEJ fixed-order generator)
* with those details
*
*/
class HepMCWriter: public EventWriter{
public:
HepMCWriter(std::string const & file, Config const & conf, LHEF::HEPRUP heprup);
~HepMCWriter() override = default;
void write(Event const & ev) override;
private:
struct HepMCWriterImpl;
struct HepMCWriterImplDeleter {
void operator()(HepMCWriterImpl* p);
};
std::unique_ptr<HepMCWriterImpl, HepMCWriterImplDeleter> impl_;
};
}
diff --git a/include/RHEJ/LesHouchesWriter.hh b/include/RHEJ/LesHouchesWriter.hh
index 1f0d63b..f8666f7 100644
--- a/include/RHEJ/LesHouchesWriter.hh
+++ b/include/RHEJ/LesHouchesWriter.hh
@@ -1,56 +1,57 @@
/** \file LesHouchesWriter.hh
* \brief Contains the writer for LesHouches Output
*/
#pragma once
#include <fstream>
-#include "RHEJ/config.hh"
#include "RHEJ/EventWriter.hh"
#include "LHEF/LHEF.h"
namespace RHEJ{
class Event;
+ struct Config;
/** \class LesHouchesWriter LesHouchesWriter.hh "include/RHEJ/LesHouchesWriter.hh"
* \brief EventWriter Class specifically for LesHouches Output
*
* This is an event writer specifically for Les Houches Output. This
* inherits from the EventWriter Class and contains its own
* contructors. TODO: in principle, this class should be movable but that
* somehow(?) breaks the write member function.
**/
class LesHouchesWriter : public EventWriter{
public:
LesHouchesWriter(std::string const & file, Config const & conf, LHEF::HEPRUP heprup);
LesHouchesWriter(LesHouchesWriter const & other) = delete;
LesHouchesWriter & operator=(LesHouchesWriter const & other) = delete;
// TODO: in principle, this class should be movable
// but that somehow(?) breaks the write member function
LesHouchesWriter(LesHouchesWriter && other) = delete;
LesHouchesWriter & operator=(LesHouchesWriter && other) = delete;
~LesHouchesWriter() override;
void write(Event const & ev) override;
private:
void write_init(){
writer_->init();
}
void rewrite_init();
LHEF::HEPRUP & heprup(){
return writer_->heprup;
}
LHEF::HEPEUP & hepeup(){
return writer_->hepeup;
}
+ Config const & config_;
std::fstream out_;
std::unique_ptr<LHEF::Writer> writer_;
};
}
diff --git a/include/RHEJ/MatrixElement.hh b/include/RHEJ/MatrixElement.hh
index f5eb2bc..ced2cb4 100644
--- a/include/RHEJ/MatrixElement.hh
+++ b/include/RHEJ/MatrixElement.hh
@@ -1,216 +1,216 @@
/** \file MatrixElement.hh
* \brief The header file which contains the MatrixElement Class
*
* This contains the MatrixElement Class which contains many functions
* used to calculate many different MatrixElements and their components.
*/
#pragma once
#include <functional>
#include "RHEJ/config.hh"
#include "RHEJ/utility.hh"
#include "RHEJ/HiggsCouplingSettings.hh"
#include "CLHEP/Vector/LorentzVector.h"
namespace RHEJ{
/** \class MatrixElement MatrixElement.hh "include/RHEJ/MatrixElement.hh
* \brief MatrixElement class. Functions for obtaining various ME and components.
*/
class MatrixElement{
public:
/** \brief MatrixElement Constructor
* @param alpha_s Function taking the renormalisation scale
* and returning the strong coupling constant
*/
MatrixElement(
std::function<double (double)> alpha_s,
Config const & conf
);
/**
* \brief regulated HEJ matrix element
* @param mur Value of the renormalisation scale
* @param incoming Incoming particles
* @param partons Outgoing particles
* @param check_momenta Special treatment for partons inside extremal jets
* @returns The HEJ matrix element including virtual corrections
*
* cf. eq. (22) in \ref Andersen:2011hs
* Incoming particles should be ordered by ascending z momentum.
* Outgoing particles should be ordered by ascending rapidity.
*
* \internal Relation to standard HEJ Met2: MatrixElement = Met2*shat^2/(pdfta*pdftb)
*/
double operator()(
double mur,
std::array<Sparticle, 2> const & incoming,
std::vector<Sparticle> const & outgoing,
bool check_momenta
) const;
//! HEJ tree-level matrix element
/**
* @param mur Value of the renormalisation scale
* @param incoming Incoming particles
* @param outgoing Outgoing particles
* @param check_momenta Special treatment for partons inside extremal jets
* @returns The HEJ matrix element without virtual corrections
*
* cf. eq. (22) in \ref Andersen:2011hs
* Incoming particles should be ordered by ascending z momentum.
* Outgoing particles should be ordered by ascending rapidity.
*/
double tree(
double mur,
std::array<Sparticle, 2> const & incoming,
std::vector<Sparticle> const & outgoing,
bool check_momenta
) const;
//! HEJ tree-level matrix element - parametric part
/**
* @param mur Value of the renormalisation scale
* @param outgoing Outgoing particles
* @returns The parametric part of the tree matrix element
*
* cf. eq. (22) in \ref Andersen:2011hs
*
* The tree level matrix element factorises into a parametric part
* which depends on the theory parameters (alpha_s and scale)
* and a kinematic part comprising the dependence on the particle momenta
* and colour factors. This function returns the former.
*/
double tree_param(
double mur,
std::array<Sparticle, 2> const & incoming,
std::vector<Sparticle> const & outgoing
) const;
//! HEJ tree-level matrix element - kinematic part
/**
* @param incoming Incoming particles
* @param partons Outgoing particles
* @param check_momenta Special treatment for partons inside extremal jets
* @returns The kinematic part of the tree matrix element
*
* cf. eq. (22) in \ref Andersen:2011hs
* Incoming particles should be ordered by ascending z momentum.
* Outgoing particles should be ordered by ascending rapidity.
*
* The tree level matrix element factorises into a parametric part
* which depends on the theory parameters (alpha_s and scale)
* and a kinematic part comprising the dependence on the particle momenta
* and colour factors. This function returns the latter.
*/
double tree_kin(
std::array<Sparticle, 2> const & incoming,
std::vector<Sparticle> const & outgoing,
bool check_momenta
) const;
/**
* \brief Calculates the Virtual Corrections
* @param mur Value of the renormalisation scale
* @param in Incoming particles
* @param out Outgoing particles
* @returns The Virtual Corrections of the Matrix Element
*
* Incoming particles should be ordered by ascending z momentum.
* Outgoing particles should be ordered by ascending rapidity.
*
* The all order virtual corrections to LL in the MRK limit is
* given by replacing 1/t in the scattering amplitude according to the
* lipatov ansatz.
*
* cf. second-to-last line of eq. (22) in \ref Andersen:2011hs
* note that indices are off by one, i.e. out[0].p corresponds to p_1
*/
double virtual_corrections(
double mur,
std::array<Sparticle, 2> const & in,
std::vector<Sparticle> const & out
) const;
private:
/**
* \brief cf. last line of eq. (22) in \ref Andersen:2011hs
* @param mur Value of Renormalisation Scale
* @param q_j ???
* @param lambda ???
*/
double omega0(
double alpha_s, double mur,
fastjet::PseudoJet const & q_j, double lambda
) const;
double tree_kin_jets(
std::array<Sparticle, 2> const & incoming,
std::vector<Sparticle> partons,
bool check_momenta
) const;
double tree_kin_Higgs(
std::array<Sparticle, 2> const & incoming,
std::vector<Sparticle> const & outgoing,
bool check_momenta
) const;
double tree_kin_Higgs_first(
std::array<Sparticle, 2> const & incoming,
std::vector<Sparticle> const & outgoing,
bool check_momenta
) const;
double tree_kin_Higgs_last(
std::array<Sparticle, 2> const & incoming,
std::vector<Sparticle> const & outgoing,
bool check_momenta
) const;
/**
* \brief Higgs inbetween extremal partons.
*
* Note that in the case of unordered emission, the Higgs is *always*
* treated as if in between the extremal (FKL) partons, even if its
* rapidity is outside the extremal parton rapidities
*/
double tree_kin_Higgs_between(
std::array<Sparticle, 2> const & incoming,
std::vector<Sparticle> const & outgoing,
bool check_momenta
) const;
double tree_param_partons(
double alpha_s, double mur,
std::vector<Sparticle> const & partons
) const;
std::vector<int> in_extremal_jet_indices(
std::vector<fastjet::PseudoJet> const & partons
) const;
std::vector<Sparticle> tag_extremal_jet_partons(
std::array<Sparticle, 2> const & incoming,
std::vector<Sparticle> out_partons, bool check_momenta
) const;
double MH2_forwardH(
RHEJ::ParticleID id,
CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in,
CLHEP::HepLorentzVector pH,
double t1, double t2
) const;
std::function<double (double)> alpha_s_;
- MartixElementParameters param;
+ MatrixElementParameters param_;
};
}
diff --git a/include/RHEJ/PhaseSpacePoint.hh b/include/RHEJ/PhaseSpacePoint.hh
index 4d32587..86cd0e5 100644
--- a/include/RHEJ/PhaseSpacePoint.hh
+++ b/include/RHEJ/PhaseSpacePoint.hh
@@ -1,157 +1,157 @@
/** \file PhaseSpacePoint.hh
* \brief Contains the PhaseSpacePoint Class
*/
#pragma once
#include <vector>
#include <CLHEP/Random/Randomize.h>
#include <CLHEP/Random/RanluxEngine.h>
#include "RHEJ/utility.hh"
#include "RHEJ/config.hh"
#include "RHEJ/Event.hh"
#include "RHEJ/Splitter.hh"
namespace RHEJ{
//! A point in resummation phase space
class PhaseSpacePoint{
public:
//! Default PhaseSpacePoint Constructor
PhaseSpacePoint() = default;
//! PhaseSpacePoint Constructor
/**
* @param ev Clustered Jet Event
* @param conf Configuration parameters
*/
PhaseSpacePoint(
Event const & ev,
Config const & conf
);
//! Get Weight Function
/**
* @returns Weight of Event
*/
double weight() const{
return weight_;
}
//! Get Incoming Function
/**
* @returns Incoming Particles
*/
std::array<Sparticle, 2> const & incoming() const{
return incoming_;
}
//! Get Outgoing Function
/**
* @returns Outgoing Particles
*/
std::vector<Sparticle> const & outgoing() const{
return outgoing_;
}
std::unordered_map<int, std::vector<Sparticle>> const & decays() const{
return decays_;
}
static constexpr int ng_max = 1000; // maximum number of extra gluons
static void reset_ranlux(std::string const & init_file);
static void reset_ranlux(char const * init_file);
private:
std::vector<fastjet::PseudoJet> cluster_jets(
std::vector<fastjet::PseudoJet> const & partons
) const;
bool pass_resummation_cuts(
std::vector<fastjet::PseudoJet> const & jets
) const;
bool pass_extremal_cuts(
fastjet::PseudoJet const & ext_parton,
fastjet::PseudoJet const & jet
) const;
int sample_ng(std::vector<fastjet::PseudoJet> const & Born_jets);
int sample_ng_jets(int ng, std::vector<fastjet::PseudoJet> const & Born_jets);
double probability_in_jet(
std::vector<fastjet::PseudoJet> const & Born_jets
) const;
std::vector<fastjet::PseudoJet> gen_non_jet(
int ng_non_jet,
double ptmin, double ptmax
);
void rescale_rapidities(
std::vector<fastjet::PseudoJet> & partons,
double ymin, double ymax
);
std::vector<fastjet::PseudoJet> reshuffle(
std::vector<fastjet::PseudoJet> const & Born_jets,
fastjet::PseudoJet const & q
);
bool jets_ok(
std::vector<fastjet::PseudoJet> const & Born_jets,
std::vector<fastjet::PseudoJet> const & partons
) const;
void reconstruct_incoming(std::array<Sparticle, 2> const & Born_incoming);
double phase_space_normalisation(
int num_Born_jets,
int num_res_partons
) const;
std::vector<fastjet::PseudoJet> split(
std::vector<fastjet::PseudoJet> const & jets, int ng_jets
);
std::vector<int> distribute_jet_partons(
int ng_jets, std::vector<fastjet::PseudoJet> const & jets
);
std::vector<fastjet::PseudoJet> split(
std::vector<fastjet::PseudoJet> const & jets,
std::vector<int> const & np_in_jet
);
bool split_preserved_jets(
std::vector<fastjet::PseudoJet> const & jets,
std::vector<fastjet::PseudoJet> const & jet_partons
) const;
template<class Particle>
Particle const & most_backward_FKL(
std::vector<Particle> const & partons
) const;
template<class Particle>
Particle const & most_forward_FKL(
std::vector<Particle> const & partons
) const;
template<class Particle>
Particle & most_backward_FKL(std::vector<Particle> & partons) const;
template<class Particle>
Particle & most_forward_FKL(std::vector<Particle> & partons) const;
bool extremal_ok(
std::vector<fastjet::PseudoJet> const & partons
) const;
void copy_AWZH_boson_from(Event const & event);
bool momentum_conserved(double ep) const;
bool unob_, unof_;
double weight_;
- const PhaseSpacePointParameters param;
+ PhaseSpacePointParameters param_;
std::array<Sparticle, 2> incoming_;
std::vector<Sparticle> outgoing_;
//! Particle decays in the format {outgoing index, decay products}
std::unordered_map<int, std::vector<Sparticle>> decays_;
static CLHEP::Ranlux64Engine ran_;
HejSplit splitter_;
};
}
diff --git a/include/RHEJ/config.hh b/include/RHEJ/config.hh
index 6b36c92..26d031c 100644
--- a/include/RHEJ/config.hh
+++ b/include/RHEJ/config.hh
@@ -1,194 +1,194 @@
/** \file config.hh
* \brief The file which handles the configuration file parameters
*
* Contains the JetParameters Struct and ScaleConfig Struct. Also
* contains the ScaleGenerator Class and the EventTreatment class.
* Config struct is also defined here. The configuration files
* parameters are read and then stored within this objects.
*/
#pragma once
#include <string>
#include <memory>
#include "fastjet/JetDefinition.hh"
#include "yaml-cpp/yaml.h"
#include "RHEJ/event_types.hh"
#include "RHEJ/Event.hh"
#include "RHEJ/HiggsCouplingSettings.hh"
#include "RHEJ/optional.hh"
#include "RHEJ/output_formats.hh"
#include "RHEJ/ScaleFunction.hh"
#include "RHEJ/utility.hh"
namespace RHEJ{
/**! \struct JetParameters config.hh "include/RHEJ/config.hh"
* \brief Contains Jet Definition and Minimum Pt to be a Jet.
*
* Contains the Fastjet definition of a Jet and a threshold (minimum) value for pt
*/
struct JetParameters{
fastjet::JetDefinition def; /**< Jet Definition */
double min_pt; /**< Minimum Jet Pt */
};
/**! \struct ScaleConfig config.hh "include/RHEJ/config.hh"
* \brief Sets up the possible choices for scale variation?
*
* There is a vector which contains the possible scale choices and
* also vector of doubles with the scale factors along with a max
* scale ratio.
*/
struct ScaleConfig{
std::vector<std::shared_ptr<ScaleFunction>> scales; /**< Vector of possible Scale choices */
std::vector<double> scale_factors; /**< Vector of possible Scale Factors */
double max_scale_ratio; /**< Maximum ratio for the scales */
};
/**! \class ScaleGenerator config.hh "include/RHEJ/config.hh"
* \brief A Class used to generate scales
*
* This class has a clustered event class and scale config struct
* within it.
*/
class ScaleGenerator{
public:
ScaleGenerator() = default;
explicit ScaleGenerator(ScaleConfig && settings):
settings_{std::move(settings)}
{}
Event operator()(Event ev) const;
ScaleConfig const & settings() const;
private:
ScaleConfig settings_;
};
/**! \enum EventTreatment config.hh "include/RHEJ/config.hh"
* \brief Enumeration which deals with how to treat events.
*
* The program will decide on how to treat an event based on
* the value of this enumeration.
*/
enum class EventTreatment{
reweight, /**< Reweigh the event */
keep, /**< Keep the event */
discard, /**< Discard the event */
};
/**
* \brief An ordered Map with EventType keys and mapped values EventTreatment
*
* If called upon with EventTreatMap.at(EventType) it will return the treatment which has
* been specified for that EventType.
*/
using EventTreatMap = std::map<event_type::EventType, EventTreatment>;
/**! \struct Config config.hh "include/RHEJ/config.hh"
* \brief Config Struct for user input parameters.
*
* The struct which handles the input card given by the user.
*
* \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
*/
struct Config {
ScaleGenerator scale_gen; /**< Scale */
JetParameters resummation_jets; /**< Resummation Jets */
JetParameters fixed_order_jets; /**< Fixed-Order Jets */
double min_extparton_pt; /**< Min External Parton Pt*/
double max_ext_soft_pt_fraction; /**< Min External Parton Pt Fraction */
int trials; /**< Number of resummation points to generate per FO */
//! Log Correct from running of \f$\alpha_s$\f On or Off
bool log_correction;
bool unweight; /**< Unweight On or Off */
std::vector<OutputFile> output; /**< Output Files Type */
optional<std::string> RanLux_init; /**< Ranlux File Choice*/
//! Map to decide what to do for differnt EventTypes
EventTreatMap treat;
YAML::Node analysis_parameters; /**< Analysis Parameters */
//! Settings for Effective Higgs-Gluon Coupling
HiggsCouplingSettings Higgs_coupling;
};
class Parameters{
protected:
- const Config & config;
+ const Config & config_;
public:
- explicit Parameters(const Config & conf): config(conf){}
- const Config & get_Config() const {return config;}
+ explicit Parameters(const Config & conf): config_(conf){}
+ const Config & get_Config() const {return config_;}
};
class PhaseSpacePointParameters: public Parameters {
public:
explicit PhaseSpacePointParameters(const Config & conf): Parameters(conf){}
///@{
/// @name getter functions
const JetParameters & jet_parameters() const
- {return config.resummation_jets;}
+ {return config_.resummation_jets;}
double jet_ptmin() const
{return jet_parameters().min_pt;}
const fastjet::JetDefinition & jet_def() const
{return jet_parameters().def;}
double jet_R() const
{return jet_def().R();}
double min_ext_pt() const
- {return config.min_extparton_pt;}
+ {return config_.min_extparton_pt;}
double max_ext_soft_pt_fraction() const
- {return config.max_ext_soft_pt_fraction;}
+ {return config_.max_ext_soft_pt_fraction;}
///@}
};
- class MartixElementParameters: public Parameters {
+ class MatrixElementParameters: public Parameters {
public:
- explicit MartixElementParameters(const Config & conf):
+ explicit MatrixElementParameters(const Config & conf):
Parameters(conf){}
///@{
/// @name getter functions
const JetParameters & jet_parameters() const
- {return config.resummation_jets;}
+ {return config_.resummation_jets;}
double jet_ptmin() const
{return jet_parameters().min_pt;}
const fastjet::JetDefinition & jet_def() const
{return jet_parameters().def;}
bool log_corr() const
- {return config.log_correction;}
+ {return config_.log_correction;}
const HiggsCouplingSettings & Hgg_setting() const
- {return config.Higgs_coupling;}
+ {return config_.Higgs_coupling;}
///@}
};
class EventReweighterParameters: public Parameters {
public:
explicit EventReweighterParameters(const Config & conf):
Parameters(conf){}
///@{
/// @name getter functions
const JetParameters & jet_parameters() const
- {return config.resummation_jets;}
+ {return config_.resummation_jets;}
double jet_ptmin() const
{return jet_parameters().min_pt;}
const fastjet::JetDefinition & jet_def() const
{return jet_parameters().def;}
const EventTreatMap & treat() const
- {return config.treat;}
+ {return config_.treat;}
///@}
};
} // namespace RHEJ
diff --git a/src/CombinedEventWriter.cc b/src/CombinedEventWriter.cc
index a8a34fa..883958d 100644
--- a/src/CombinedEventWriter.cc
+++ b/src/CombinedEventWriter.cc
@@ -1,22 +1,24 @@
#include "RHEJ/CombinedEventWriter.hh"
+#include "RHEJ/config.hh"
+
namespace RHEJ{
CombinedEventWriter::CombinedEventWriter(
Config const & conf,
LHEF::HEPRUP const & heprup
- ):EventWriter(conf){
+ ){
std::vector<OutputFile> const & outfiles = conf.output;
writers_.reserve(outfiles.size());
for(OutputFile const & outfile: outfiles){
writers_.emplace_back(
make_format_writer(outfile.format, outfile.name, conf, heprup)
);
}
}
void CombinedEventWriter::write(Event const & ev){
for(auto & writer: writers_) writer->write(ev);
}
}
diff --git a/src/EventReweighter.cc b/src/EventReweighter.cc
index 3bdf70e..4645ab4 100644
--- a/src/EventReweighter.cc
+++ b/src/EventReweighter.cc
@@ -1,302 +1,301 @@
#include "RHEJ/EventReweighter.hh"
#include <string>
#include <unordered_map>
#include "RHEJ/PhaseSpacePoint.hh"
#include "RHEJ/PDG_codes.hh"
-#include "RHEJ/config.hh" // ScaleGenerator definition
#include "RHEJ/debug.hh"
namespace RHEJ{
using EventType = event_type::EventType;
namespace {
static_assert(
std::numeric_limits<double>::has_quiet_NaN,
"no quiet NaN for double"
);
constexpr double NaN = std::numeric_limits<double>::quiet_NaN();
UnclusteredEvent to_UnclusteredEvent(PhaseSpacePoint const & psp){
UnclusteredEvent result;
result.incoming = psp.incoming();
std::sort(
begin(result.incoming), end(result.incoming),
[](Sparticle o1, Sparticle o2){return o1.p.pz()<o2.p.pz();}
);
assert(result.incoming.size() == 2);
result.outgoing = psp.outgoing();
assert(
std::is_sorted(
begin(result.outgoing), end(result.outgoing),
rapidity_less{}
)
);
assert(result.outgoing.size() >= 2);
result.decays = psp.decays();
result.central.mur = NaN;
result.central.muf = NaN;
result.central.weight = psp.weight();
return result;
}
} // namespace anonymous
EventReweighter::EventReweighter(
LHEF::HEPRUP const & heprup,
Config const & conf
):
EventReweighter{
RHEJ::Beam{
heprup.EBMUP.first,
{{
static_cast<RHEJ::ParticleID>(heprup.IDBMUP.first),
static_cast<RHEJ::ParticleID>(heprup.IDBMUP.second)
}}
},
heprup.PDFSUP.first,
conf
}
{
if(heprup.EBMUP.second != E_beam_){
throw std::invalid_argument(
"asymmetric beam: " + std::to_string(E_beam_)
+ " ---> <--- " + std::to_string(heprup.EBMUP.second)
);
};
if(heprup.PDFSUP.second != pdf_.id()){
throw std::invalid_argument(
"conflicting PDF ids: " + std::to_string(pdf_.id())
+ " vs. " + std::to_string(heprup.PDFSUP.second)
);
}
}
EventReweighter::EventReweighter(
Beam beam,
int pdf_id,
Config const & conf
):
- param(conf),
+ param_(conf),
E_beam_{beam.E},
pdf_{pdf_id, beam.type.front(), beam.type.back()},
MEt2_{
[this](double mu){ return pdf_.Halphas(mu); },
- param.get_Config()
+ param_.get_Config()
}
{}
PDF const & EventReweighter::pdf() const{
return pdf_;
}
namespace{
static_assert(
std::numeric_limits<double>::has_infinity, "infinity not supported"
);
// there is no simple way to create a vector of move-only objects
// this is a clumsy work-around
ScaleGenerator make_identity(){
std::vector<std::shared_ptr<ScaleFunction>> scale_fun;
scale_fun.emplace_back(new InputScales);
return ScaleGenerator{
ScaleConfig{
std::move(scale_fun), {}, std::numeric_limits<double>::infinity()
}
};
}
const ScaleGenerator identity{make_identity()};
}
std::vector<Event> EventReweighter::reweight(
Event const & input_ev, int num_events
){
return reweight(input_ev, num_events, identity);
}
std::vector<Event> EventReweighter::reweight(
Event const & input_ev, int num_events,
ScaleGenerator const & gen_scales
){
auto res_events = gen_res_events(input_ev, num_events);
if(res_events.empty()) return {};
for(auto & event: res_events) event = gen_scales(event);
return rescale(input_ev, std::move(res_events));
}
/**
* \brief main generation/reweighting function: generate phase space points and divide out Born factors
*/
std::vector<Event> EventReweighter::gen_res_events(
Event const & ev,
int phase_space_points
){
assert(ev.variations().empty());
- switch(param.treat().at(ev.type())){
+ switch(param_.treat().at(ev.type())){
case EventTreatment::discard: return {};
case EventTreatment::keep:
if(! jets_pass_resummation_cuts(ev)) return {};
else return {ev};
default:;
}
const double Born_shat = shat(ev);
std::vector<Event> resummation_events;
for(int psp_number = 0; psp_number < phase_space_points; ++psp_number){
PhaseSpacePoint psp{
ev,
- param.get_Config()
+ param_.get_Config()
};
if(psp.weight() == 0.) continue;
if(psp.incoming()[0].E() > E_beam_ || psp.incoming()[1].E() > E_beam_) continue;
resummation_events.emplace_back(
- to_UnclusteredEvent(std::move(psp)), param.jet_def(), param.jet_ptmin()
+ to_UnclusteredEvent(std::move(psp)), param_.jet_def(), param_.jet_ptmin()
);
auto & new_event = resummation_events.back();
assert(new_event.variations().empty());
new_event.central().mur = ev.central().mur;
new_event.central().muf = ev.central().muf;
const double resum_shat = shat(new_event);
new_event.central().weight *= ev.central().weight*Born_shat*Born_shat/
(phase_space_points*resum_shat*resum_shat);
}
return resummation_events;
}
std::vector<Event> EventReweighter::rescale(
Event const & Born_ev,
std::vector<Event> events
) const{
const double Born_pdf = pdf_factors(Born_ev).central;
const double Born_ME = tree_matrix_element(Born_ev);
for(auto & cur_event: events){
const auto pdf = pdf_factors(cur_event);
assert(pdf.variations.size() == cur_event.variations().size());
const auto ME = matrix_elements(cur_event);
assert(ME.variations.size() == cur_event.variations().size());
cur_event.central().weight *= pdf.central*ME.central/(Born_pdf*Born_ME);
for(size_t i = 0; i < cur_event.variations().size(); ++i){
cur_event.variations(i).weight *=
pdf.variations[i]*ME.variations[i]/(Born_pdf*Born_ME);
}
}
return events;
};
bool EventReweighter::jets_pass_resummation_cuts(
Event const & ev
) const{
const auto out_as_PseudoJet = to_PseudoJet(filter_partons(ev.outgoing()));
- fastjet::ClusterSequence cs{out_as_PseudoJet, param.jet_def()};
- return cs.inclusive_jets(param.jet_ptmin()).size() == ev.jets().size();
+ fastjet::ClusterSequence cs{out_as_PseudoJet, param_.jet_def()};
+ return cs.inclusive_jets(param_.jet_ptmin()).size() == ev.jets().size();
}
EventReweighter::EventFactors
EventReweighter::pdf_factors(Event const & ev) const{
auto const & a = ev.incoming().front();
auto const & b = ev.incoming().back();
const double xa = a.p.e()/E_beam_;
const double xb = b.p.e()/E_beam_;
EventFactors result;
std::unordered_map<double, double> known_pdf;
result.central =
pdf_.pdfpt(0,xa,ev.central().muf,a.type)*
pdf_.pdfpt(1,xb,ev.central().muf,b.type);
known_pdf.emplace(ev.central().muf, result.central);
result.variations.reserve(ev.variations().size());
for(auto const & ev_param: ev.variations()){
const double muf = ev_param.muf;
auto cur_pdf = known_pdf.find(muf);
if(cur_pdf == known_pdf.end()){
cur_pdf = known_pdf.emplace(
muf,
pdf_.pdfpt(0,xa,muf,a.type)*pdf_.pdfpt(1,xb,muf,b.type)
).first;
}
result.variations.emplace_back(cur_pdf->second);
}
assert(result.variations.size() == ev.variations().size());
return result;
}
EventReweighter::EventFactors
EventReweighter::matrix_elements(Event const & ev) const{
- assert(param.treat().count(ev.type()) > 0);
- if(param.treat().find(ev.type())->second == EventTreatment::keep){
+ assert(param_.treat().count(ev.type()) > 0);
+ if(param_.treat().find(ev.type())->second == EventTreatment::keep){
return fixed_order_scale_ME(ev);
}
// precompute overall kinematic factor
const double ME_kin = MEt2_.tree_kin(ev.incoming(), ev.outgoing(), true);
EventFactors result;
std::unordered_map<double, double> known_ME;
result.central = MEt2_(
ev.central().mur,
ev.incoming(), ev.outgoing(),
true
);
known_ME.emplace(ev.central().mur, result.central);
result.variations.reserve(ev.variations().size());
- for(auto const & param: ev.variations()){
- const double mur = param.mur;
+ for(auto const & param_: ev.variations()){
+ const double mur = param_.mur;
auto cur_ME = known_ME.find(mur);
if(cur_ME == known_ME.end()){
const double ME = MEt2_.tree_param(
mur, ev.incoming(), ev.outgoing()
)*ME_kin*MEt2_.virtual_corrections(
mur, ev.incoming(), ev.outgoing()
);
cur_ME = known_ME.emplace(mur, ME).first;
}
result.variations.emplace_back(cur_ME->second);
}
assert(result.variations.size() == ev.variations().size());
return result;
}
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){
+ assert(param_.treat().count(ev.type()) > 0);
+ if(param_.treat().find(ev.type())->second == EventTreatment::keep){
return fixed_order_scale_ME(ev).central;
}
return MEt2_.tree(
ev.central().mur,
ev.incoming(), ev.outgoing(),
false
);
}
EventReweighter::EventFactors
EventReweighter::fixed_order_scale_ME(Event const & ev) const{
const int alpha_s_power = std::count_if(
begin(ev.outgoing()), end(ev.outgoing()),
[](Sparticle const & p){ return is_parton(p); }
);
EventFactors result;
result.central = pow(pdf_.Halphas(ev.central().mur), alpha_s_power);
for(auto const & var: ev.variations()){
result.variations.emplace_back(
pow(pdf_.Halphas(var.mur), alpha_s_power)
);
}
return result;
}
}
diff --git a/src/HepMCWriter.cc b/src/HepMCWriter.cc
index b1eb236..f9820d8 100644
--- a/src/HepMCWriter.cc
+++ b/src/HepMCWriter.cc
@@ -1,242 +1,243 @@
#include "RHEJ/HepMCWriter.hh"
#include <cassert>
#ifdef RHEJ_BUILD_WITH_HepMC_VERSION
#if RHEJ_BUILD_WITH_HepMC_VERSION >= 3
#include "HepMC/WriterAscii.h"
#include "HepMC/LHEFAttributes.h"
#else
#include "HepMC/IO_GenEvent.h"
#endif
#include "HepMC/GenVertex.h"
#include "HepMC/GenParticle.h"
+#include "RHEJ/config.hh"
#include "RHEJ/utility.hh"
#include "RHEJ/Event.hh"
#include "RHEJ/exceptions.hh"
#include "RHEJ/Version.hh"
namespace RHEJ{
#if RHEJ_BUILD_WITH_HepMC_VERSION >= 3
namespace {
void add_generator_tag(LHEF::HEPRUP & heprup){
heprup.generators.emplace_back(LHEF::XMLTag{});
heprup.generators.back().name = package_name();
heprup.generators.back().version = version();
}
void reset_weight_info(LHEF::HEPRUP & heprup){
heprup.IDWTUP = -4;
// use placeholders for unknown init block values
// we can overwrite them after processing all events
heprup.XSECUP = {0.};
heprup.XERRUP = {0.};
heprup.XMAXUP = {0.};
}
HepMC::FourVector to_FourVector(Sparticle const & sp){
return {sp.px(), sp.py(), sp.pz(), sp.E()};
}
HepMC::shared_ptr<HepMC::GenRunInfo> & add_HEPRUP(
HepMC::shared_ptr<HepMC::GenRunInfo> & runinfo,
LHEF::HEPRUP heprup
){
add_generator_tag(heprup);
reset_weight_info(heprup);
auto hepr = HepMC::make_shared<HepMC::HEPRUPAttribute>();
hepr->heprup = std::move(heprup);
runinfo->add_attribute("HEPRUP", hepr);
for (int i = 0, N = hepr->heprup.generators.size(); i < N; ++i ){
HepMC::GenRunInfo::ToolInfo tool;
tool.name = hepr->heprup.generators[i].name;
tool.version = hepr->heprup.generators[i].version;
tool.description = hepr->heprup.generators[i].contents;
runinfo->tools().push_back(tool);
}
return runinfo;
}
} // namespace anonymous
struct HepMCWriter::HepMCWriterImpl: public EventWriter{
+ Config config_;
HepMC::shared_ptr<HepMC::GenRunInfo> runinfo_;
HepMC::WriterAscii writer_;
HepMCWriterImpl(
std::string const & file, Config const & conf, LHEF::HEPRUP && heprup
):
- EventWriter(conf),
+ config_(conf),
runinfo_{HepMC::make_shared<HepMC::GenRunInfo>()},
writer_{file, add_HEPRUP(runinfo_, std::move(heprup))}
{}
HepMCWriterImpl & operator=(HepMCWriterImpl const & other) = delete;
HepMCWriterImpl(HepMCWriterImpl const & other) = delete;
HepMCWriterImpl & operator=(HepMCWriterImpl && other) = delete;
HepMCWriterImpl(HepMCWriterImpl && other) = delete;
~HepMCWriterImpl(){
writer_.close();
}
void write(Event const & ev){
if(!ev.decays().empty()){
throw not_implemented{"HepMC output for events with decays"};
}
auto vx = HepMC::make_shared<HepMC::GenVertex>();
for(auto const & in: ev.incoming()){
vx->add_particle_in(
HepMC::make_shared<HepMC::GenParticle>(
to_FourVector(in), static_cast<int>(in.type), -1
)
);
}
for(auto const & out: ev.outgoing()){
vx->add_particle_out(
HepMC::make_shared<HepMC::GenParticle>(
to_FourVector(out), static_cast<int>(out.type), +1
)
);
}
HepMC::GenEvent out_ev{runinfo_, HepMC::Units::GEV, HepMC::Units::MM};
out_ev.add_vertex(vx);
out_ev.weights().reserve(ev.variations().size() + 1u);
out_ev.weights().emplace_back(ev.central().weight);
for(auto const & var: ev.variations()){
out_ev.weights().emplace_back(var.weight);
}
writer_.write_event(out_ev);
}
};
void HepMCWriter::HepMCWriterImplDeleter::operator()(HepMCWriterImpl* p) {
delete p;
}
#else // HepMC 2
namespace{
HepMC::FourVector to_FourVector(Sparticle const & sp){
return {sp.px(), sp.py(), sp.pz(), sp.E()};
}
}
struct HepMCWriter::HepMCWriterImpl: public EventWriter{
struct HepMCInfo{
size_t event_count = 0.;
double tot_weight = 0.;
double tot_weight2 = 0.;
void add_event(Event const & ev){
double wt = ev.central().weight;
tot_weight += wt;
tot_weight2 += wt * wt;
++event_count;
}
};
+ Config config_;
HepMC::IO_GenEvent writer_;
HepMCInfo genInfo_;
HepMCWriterImpl(
std::string const & file, Config const & conf, LHEF::HEPRUP &&
):
- EventWriter(conf),
+ config_(conf),
writer_{file}
{}
HepMCWriterImpl & operator=(HepMCWriterImpl const & other) = delete;
HepMCWriterImpl(HepMCWriterImpl const & other) = delete;
HepMCWriterImpl & operator=(HepMCWriterImpl && other) = delete;
HepMCWriterImpl(HepMCWriterImpl && other) = delete;
void write(Event const & ev){
if(!ev.decays().empty()){
throw not_implemented{"HepMC output for events with decays"};
}
auto vx = new HepMC::GenVertex();
for(auto const & in: ev.incoming()){
vx->add_particle_in(
new HepMC::GenParticle{
to_FourVector(in), static_cast<int>(in.type), -1
}
);
}
for(auto const & out: ev.outgoing()){
vx->add_particle_out(
new HepMC::GenParticle{
to_FourVector(out), static_cast<int>(out.type), +1
}
);
}
HepMC::GenEvent out_ev{HepMC::Units::GEV, HepMC::Units::MM};
out_ev.add_vertex(vx);
/// weights
out_ev.weights().push_back(ev.central().weight);
for(auto const & var: ev.variations()){
out_ev.weights().push_back(var.weight);
}
/// @TODO add name list
/// general informations
genInfo_.add_event(ev);
out_ev.set_signal_process_id(ev.type()+1); /// "+1": conistent with lhe
out_ev.set_event_scale(ev.central().mur); /// event scale
out_ev.set_event_number(genInfo_.event_count); /// event number
/// cross section
HepMC::GenCrossSection t_xs;
t_xs.set_cross_section( genInfo_.tot_weight , sqrt(genInfo_.tot_weight2) );
out_ev.set_cross_section( t_xs );
/// @TODO add alphaQCD (need function) and alphaQED
/// @TODO output pdf (currently not avaiable from event alone)
writer_.write_event(&out_ev);
}
};
void HepMCWriter::HepMCWriterImplDeleter::operator()(HepMCWriterImpl* p) {
delete p;
}
#endif // HepMC 2
HepMCWriter::HepMCWriter(std::string const & file, Config const & conf, LHEF::HEPRUP heprup):
- EventWriter(conf),
impl_{std::unique_ptr<HepMCWriterImpl, HepMCWriterImplDeleter>{
new HepMCWriterImpl(file, conf, std::move(heprup))
}}
{}
void HepMCWriter::write(Event const & ev){
impl_->write(ev);
}
} // namespace RHEJ
#else // no HepMC
namespace RHEJ{
class HepMCWriter::HepMCWriterImpl{};
- HepMCWriter::HepMCWriter(std::string const &, Config const & conf, LHEF::HEPRUP):
- EventWriter(conf){
+ HepMCWriter::HepMCWriter(std::string const &, Config const & conf, LHEF::HEPRUP){
throw std::invalid_argument(
"Failed to create HepMC writer: "
"Reversed HEJ was built without HepMC support"
);
}
void HepMCWriter::write(Event const &){
assert(false);
}
void HepMCWriter::HepMCWriterImplDeleter::operator()(HepMCWriterImpl* p) {
delete p;
}
}
#endif
diff --git a/src/LesHouchesWriter.cc b/src/LesHouchesWriter.cc
index 8532b41..ad85e97 100644
--- a/src/LesHouchesWriter.cc
+++ b/src/LesHouchesWriter.cc
@@ -1,91 +1,92 @@
#include <stdexcept>
#include <memory>
#include <cassert>
#include "RHEJ/LesHouchesWriter.hh"
#include "RHEJ/event_types.hh"
+#include "RHEJ/config.hh"
#include "RHEJ/Event.hh"
#include "RHEJ/Version.hh"
namespace RHEJ{
namespace{
template<class T, class... Args>
std::unique_ptr<T> make_unique(Args&&... a){
return std::unique_ptr<T>{new T{std::forward<Args>(a)...}};
}
}
LesHouchesWriter::LesHouchesWriter(
std::string const & file, Config const & conf, LHEF::HEPRUP heprup
):
- EventWriter(conf),
+ config_(conf),
out_{file, std::fstream::in | std::fstream::out | std::fstream::trunc},
writer_{RHEJ::make_unique<LHEF::Writer>(out_)}
{
if(! out_.is_open()){
throw std::ios_base::failure("Failed to open " + file);
};
writer_->heprup = std::move(heprup);
// lhe Stardard: IDWTUP (negative => weights = +/-)
// 3: weight=+/-, xs given in head (same as default MG)
// 4: weight=+/-, xs = avg(weights)
writer_->heprup.IDWTUP = -3;
writer_->heprup.generators.emplace_back(LHEF::XMLTag{});
writer_->heprup.generators.back().name = package_name();
writer_->heprup.generators.back().version = version();
// use placeholders for unknown init block values
// we can overwrite them after processing all events
writer_->heprup.XSECUP = std::vector<double>(event_type::last_type+1, 0.);
writer_->heprup.XERRUP = std::vector<double>(event_type::last_type+1, 0.);
writer_->heprup.XMAXUP = std::vector<double>(event_type::last_type+1, 0.);
write_init();
}
void LesHouchesWriter::write(Event const & ev){
assert(writer_ && out_.is_open());
const double wt = ev.central().weight;
writer_->hepeup = RHEJ::to_HEPEUP(std::move(ev), &heprup());
writer_->writeEvent();
heprup().XSECUP[ev.type()] += wt;
heprup().XERRUP[ev.type()] += wt*wt;
if(wt > heprup().XMAXUP[ev.type()]){
heprup().XMAXUP[ev.type()] = wt;
}
}
// this function is called after overwritting the Les Houches init block
// assert that we have overwritten *exactly* the init block,
// i.e. an intact event block is next
void assert_next_event_intact(std::istream & out){
(void) out; // suppress compiler warnings if not in debug mode
#ifndef NDEBUG
std::string line;
getline(out, line);
assert(line == "<event>");
#endif
}
void LesHouchesWriter::rewrite_init(){
assert(writer_ && out_.is_open());
// replace placeholder entries
const auto pos = out_.tellp();
out_.seekp(0);
writer_->init();
assert_next_event_intact(out_);
out_.seekp(pos);
}
LesHouchesWriter::~LesHouchesWriter(){
assert(writer_ && out_.is_open());
for(auto & xs_err: heprup().XERRUP)
{
xs_err = sqrt(xs_err);
}
rewrite_init();
}
}
diff --git a/src/MatrixElement.cc b/src/MatrixElement.cc
index dd636dd..8f89451 100644
--- a/src/MatrixElement.cc
+++ b/src/MatrixElement.cc
@@ -1,748 +1,748 @@
#include "RHEJ/MatrixElement.hh"
#include <CLHEP/Random/Randomize.h>
#include <CLHEP/Random/RanluxEngine.h>
#include "RHEJ/Constants.hh"
#include "RHEJ/currents.hh"
#include "RHEJ/PDG_codes.hh"
#include "RHEJ/uno.hh"
#include "RHEJ/debug.hh"
namespace RHEJ{
//cf. last line of eq. (22) in \ref Andersen:2011hs
double MatrixElement::omega0(
double alpha_s, double mur,
fastjet::PseudoJet const & q_j, double lambda
) const {
const double result = - alpha_s*N_C/M_PI*log(q_j.perp2()/(lambda*lambda));
- if(! param.log_corr()) return result;
+ if(! param_.log_corr()) return result;
// use alpha_s(sqrt(q_j*lambda)), evolved to mur
return (
1. + alpha_s/(4.*M_PI)*beta0*log(mur*mur/(q_j.perp()*lambda))
)*result;
}
double MatrixElement::virtual_corrections(
double mur,
std::array<Sparticle, 2> const & in,
std::vector<Sparticle> const & out
) const{
fastjet::PseudoJet const & pa = in.front().p;
#ifndef NDEBUG
fastjet::PseudoJet const & pb = in.back().p;
#endif
assert(std::is_sorted(out.begin(), out.end(), rapidity_less{}));
assert(out.size() >= 2);
assert(pa.pz() < pb.pz());
fastjet::PseudoJet q = pa - out[0].p;
size_t first_idx = 0;
size_t last_idx = out.size() - 1;
// if there is a Higgs or unordered gluon outside the extremal partons
// then it is not part of the FKL ladder and does not contribute
// to the virtual corrections
if(out.front().type == pid::Higgs || has_unob_gluon(in, out)){
q -= out[1].p;
++first_idx;
}
if(out.back().type == pid::Higgs || has_unof_gluon(in, out)){
--last_idx;
}
double exponent = 0;
const double alpha_s = alpha_s_(mur);
for(size_t j = first_idx; j < last_idx; ++j){
exponent += omega0(alpha_s, mur, q, CLAMBDA)*(
out[j+1].rapidity() - out[j].rapidity()
);
q -= out[j+1].p;
}
assert(
nearby(q, -1*pb)
|| out.back().type == pid::Higgs
|| has_unof_gluon(in, out)
);
return exp(exponent);
}
}
namespace {
//! Lipatov vertex for partons emitted into extremal jets
double C2Lipatov(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2)
{
CLHEP::HepLorentzVector temptrans=-(qav+qbv);
CLHEP::HepLorentzVector p5=qav-qbv;
CLHEP::HepLorentzVector CL=temptrans+(qav.m2()/p5.dot(p1)+2.*p5.dot(p2)/p1.dot(p2))*p1-p2*(qbv.m2()/p5.dot(p2)+2.*p5.dot(p1)/p1.dot(p2));
#if printoutput
cout << "#Fadin qa : "<<qav<<endl;
cout << "#Fadin qb : "<<qbv<<endl;
cout << "#Fadin p1 : "<<p1<<endl;
cout << "#Fadin p2 : "<<p2<<endl;
cout << "#Fadin p5 : "<<p5<<endl;
cout << "#Fadin Gauge Check : "<< CL.dot(p5)<<endl;
cout << "#Fadin C2L : "<< -CL.dot(CL)<<" "<<-CL.dot(CL)/(qav.m2()*qbv.m2())/(4./p5.perp2())<<endl;
#endif
#if 0
if (-CL.dot(CL)<0.)
// if (fabs(CL.dot(p5))>fabs(CL.dot(CL))) // not sufficient!
return 0.;
else
#endif
return -CL.dot(CL);
}
//! Lipatov vertex with soft subtraction for partons emitted into extremal jets
double C2Lipatovots(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2)
{
double kperp=(qav-qbv).perp();
if (kperp>RHEJ::CLAMBDA)
return C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2());
else {
double Cls=(C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2()));
double temp=Cls-4./(kperp*kperp);
// std::cout <<kperp <<" "<<temp<<" "<<4./(kperp*kperp)<<" "<<(C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2()))<<std::endl;
return temp;
}
}
//! Lipatov vertex
double C2Lipatov(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv, CLHEP::HepLorentzVector pim, CLHEP::HepLorentzVector pip,CLHEP::HepLorentzVector pom, CLHEP::HepLorentzVector pop) // B
{
CLHEP::HepLorentzVector temptrans=-(qav+qbv);
CLHEP::HepLorentzVector p5=qav-qbv;
CLHEP::HepLorentzVector CL=temptrans+qav.m2()*(1./p5.dot(pip)*pip+1./p5.dot(pop)*pop)/2.-qbv.m2()*(1./p5.dot(pim)*pim+1./p5.dot(pom)*pom)/2.+(pip*(p5.dot(pim)/pip.dot(pim)+p5.dot(pom)/pip.dot(pom)) + pop*(p5.dot(pim)/pop.dot(pim)+p5.dot(pom)/pop.dot(pom)) - pim*(p5.dot(pip)/pip.dot(pim) + p5.dot(pop)/pop.dot(pim)) - pom*(p5.dot(pip)/pip.dot(pom) + p5.dot(pop)/pop.dot(pom)) )/2.;
return -CL.dot(CL);
}
//! Lipatov vertex with soft subtraction
double C2Lipatovots(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv, CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pb, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2)
{
double kperp=(qav-qbv).perp();
if (kperp>RHEJ::CLAMBDA)
return C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2());
else {
double Cls=(C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2()));
double temp=Cls-4./(kperp*kperp);
return temp;
}
}
/** Matrix element squared for tree-level current-current scattering
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @returns ME Squared for Tree-Level Current-Current Scattering
*/
double ME_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa
){
if (aptype==21&&bptype==21) {
return jM2gg(pn,pb,p1,pa);
} else if (aptype==21&&bptype!=21) {
if (bptype > 0)
return jM2qg(pn,pb,p1,pa);
else
return jM2qbarg(pn,pb,p1,pa);
}
else if (bptype==21&&aptype!=21) { // ----- || -----
if (aptype > 0)
return jM2qg(p1,pa,pn,pb);
else
return jM2qbarg(p1,pa,pn,pb);
}
else { // they are both quark
if (bptype>0) {
if (aptype>0)
return jM2qQ(pn,pb,p1,pa);
else
return jM2qQbar(pn,pb,p1,pa);
}
else {
if (aptype>0)
return jM2qQbar(p1,pa,pn,pb);
else
return jM2qbarQbar(pn,pb,p1,pa);
}
}
throw std::logic_error("unknown particle types");
}
/** \brief Matrix element squared for tree-level current-current scattering with Higgs
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param qH t-channel momentum before Higgs
* @param qHp1 t-channel momentum after Higgs
* @returns ME Squared for Tree-Level Current-Current Scattering with Higgs
*/
double ME_Higgs_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
double mt, bool include_bottom, double mb
){
if (aptype==21&&bptype==21) // gg initial state
return MH2gg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else if (aptype==21&&bptype!=21) {
if (bptype > 0)
return MH2qg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4./9.;
else
return MH2qbarg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4./9.;
}
else if (bptype==21&&aptype!=21) {
if (aptype > 0)
return MH2qg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4./9.;
else
return MH2qbarg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4./9.;
}
else { // they are both quark
if (bptype>0) {
if (aptype>0)
return MH2qQ(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.);
else
return MH2qQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.);
}
else {
if (aptype>0)
return MH2qQbar(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4.*4./(9.*9.);
else
return MH2qbarQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.);
}
}
throw std::logic_error("unknown particle types");
}
/** \brief Current matrix element squared with Higgs and unordered forward emission
* @param aptype Particle A PDG ID
* @param bptype Particle B PDG ID
* @param punof Unordered Particle Momentum
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param qH t-channel momentum before Higgs
* @param qHp1 t-channel momentum after Higgs
* @returns ME Squared with Higgs and unordered forward emission
*/
double ME_Higgs_current_unof(
int aptype, int bptype,
CLHEP::HepLorentzVector const & punof,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
double mt, bool include_bottom, double mb
){
if (aptype==21&&bptype!=21) {
if (bptype > 0)
return jM2unogqHg(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unogqbarHg(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
else { // they are both quark
if (bptype>0) {
if (aptype>0)
return jM2unogqHQ(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unogqHQbar(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
else {
if (aptype>0)
return jM2unogqbarHQ(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unogqbarHQbar(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
}
throw std::logic_error("unknown particle types");
}
/** \brief Current matrix element squared with Higgs and unordered backward emission
* @param aptype Particle A PDG ID
* @param bptype Particle B PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param punob Unordered back Particle Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param qH t-channel momentum before Higgs
* @param qHp1 t-channel momentum after Higgs
* @returns ME Squared with Higgs and unordered backward emission
*/
double ME_Higgs_current_unob(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & punob,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
double mt, bool include_bottom, double mb
){
if (bptype==21&&aptype!=21) {
if (aptype > 0)
return jM2unobgHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unobgHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
else { // they are both quark
if (aptype>0) {
if (bptype>0)
return jM2unobqHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unobqbarHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
else {
if (bptype>0)
return jM2unobqHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unobqbarHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
}
throw std::logic_error("unknown particle types");
}
CLHEP::HepLorentzVector to_HepLorentzVector(RHEJ::Sparticle const & particle){
return {particle.p.px(), particle.p.py(), particle.p.pz(), particle.p.E()};
}
} // namespace RHEJ
namespace RHEJ{
MatrixElement::MatrixElement(
std::function<double (double)> alpha_s,
Config const & conf
):
alpha_s_{std::move(alpha_s)},
- param(conf)
+ param_(conf)
{}
double MatrixElement::operator()(
double mur,
std::array<Sparticle, 2> const & incoming,
std::vector<Sparticle> const & outgoing,
bool check_momenta
) const {
return tree(
mur,
incoming, outgoing,
check_momenta
)*virtual_corrections(
mur,
incoming, outgoing
);
}
double MatrixElement::tree_kin(
std::array<Sparticle, 2> const & incoming,
std::vector<Sparticle> const & outgoing,
bool check_momenta
) const {
assert(
std::is_sorted(
incoming.begin(), incoming.end(),
[](Sparticle o1, Sparticle o2){return o1.p.pz()<o2.p.pz();}
)
);
assert(std::is_sorted(outgoing.begin(), outgoing.end(), rapidity_less{}));
auto AWZH_boson = std::find_if(
begin(outgoing), end(outgoing),
[](Sparticle const & p){return is_AWZH_boson(p);}
);
if(AWZH_boson == end(outgoing)){
return tree_kin_jets(incoming, outgoing, check_momenta);
}
switch(AWZH_boson->type){
case pid::Higgs: {
static constexpr double mH = 125.;
const double alpha_s_mH = alpha_s_(mH);
return alpha_s_mH*alpha_s_mH/(256.*pow(M_PI, 5))*tree_kin_Higgs(
incoming, outgoing, check_momenta
);
}
// TODO
case pid::Wp:
case pid::Wm:
case pid::photon:
case pid::Z:
default:
throw std::logic_error("Emission of boson of unsupported type");
}
}
namespace{
constexpr int extremal_jet_idx = 1;
constexpr int no_extremal_jet_idx = 0;
bool treat_as_extremal(Sparticle const & parton){
return parton.p.user_index() == extremal_jet_idx;
}
template<class InputIterator>
double FKL_ladder_weight(
InputIterator begin_gluon, InputIterator end_gluon,
CLHEP::HepLorentzVector const & q0,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pn
){
double wt = 1;
auto qi = q0;
for(auto gluon_it = begin_gluon; gluon_it != end_gluon; ++gluon_it){
assert(gluon_it->type == pid::gluon);
const auto g = to_HepLorentzVector(*gluon_it);
const auto qip1 = qi - g;
if(treat_as_extremal(*gluon_it)){
wt *= C2Lipatovots(qip1, qi, pa, pb)*C_A;
} else{
wt *= C2Lipatovots(qip1, qi, pa, pb, p1, pn)*C_A;
}
qi = qip1;
}
return wt;
}
} // namespace anonymous
std::vector<Sparticle> MatrixElement::tag_extremal_jet_partons(
std::array<Sparticle, 2> const & incoming,
std::vector<Sparticle> out_partons, bool check_momenta
) const{
if(!check_momenta){
for(auto & parton: out_partons){
parton.p.set_user_index(no_extremal_jet_idx);
}
return out_partons;
}
- fastjet::ClusterSequence cs(to_PseudoJet(out_partons), param.jet_def());
- const auto jets = sorted_by_rapidity(cs.inclusive_jets(param.jet_ptmin()));
+ fastjet::ClusterSequence cs(to_PseudoJet(out_partons), param_.jet_def());
+ const auto jets = sorted_by_rapidity(cs.inclusive_jets(param_.jet_ptmin()));
assert(jets.size() >= 2);
auto most_backward = begin(jets);
auto most_forward = end(jets) - 1;
// skip jets caused by unordered emission
if(has_unob_gluon(incoming, out_partons)){
assert(jets.size() >= 3);
++most_backward;
}
else if(has_unof_gluon(incoming, out_partons)){
assert(jets.size() >= 3);
--most_forward;
}
const auto extremal_jet_indices = cs.particle_jet_indices(
{*most_backward, *most_forward}
);
assert(extremal_jet_indices.size() == out_partons.size());
for(size_t i = 0; i < out_partons.size(); ++i){
assert(RHEJ::is_parton(out_partons[i]));
const int idx = (extremal_jet_indices[i]>=0)?
extremal_jet_idx:
no_extremal_jet_idx;
out_partons[i].p.set_user_index(idx);
}
return out_partons;
}
double MatrixElement::tree_kin_jets(
std::array<Sparticle, 2> const & incoming,
std::vector<Sparticle> partons,
bool check_momenta
) const {
partons = tag_extremal_jet_partons(incoming, partons, check_momenta);
if(has_unob_gluon(incoming, partons) || has_unof_gluon(incoming, partons)){
throw std::logic_error("unordered emission not implemented for pure jets");
}
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto p1 = to_HepLorentzVector(partons.front());
const auto pn = to_HepLorentzVector(partons.back());
return ME_current(
incoming[0].type, incoming[1].type,
pn, pb, p1, pa
)/(4*(N_C*N_C - 1))*FKL_ladder_weight(
begin(partons) + 1, end(partons) - 1,
pa - p1, pa, pb, p1, pn
);
}
double MatrixElement::tree_kin_Higgs(
std::array<Sparticle, 2> const & incoming,
std::vector<Sparticle> const & outgoing,
bool check_momenta
) const {
if(has_uno_gluon(incoming, outgoing)){
return tree_kin_Higgs_between(incoming, outgoing, check_momenta);
}
if(outgoing.front().type == pid::Higgs){
return tree_kin_Higgs_first(incoming, outgoing, check_momenta);
}
if(outgoing.back().type == pid::Higgs){
return tree_kin_Higgs_last(incoming, outgoing, check_momenta);
}
return tree_kin_Higgs_between(incoming, outgoing, check_momenta);
}
double MatrixElement::MH2_forwardH(
RHEJ::ParticleID id,
CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in,
CLHEP::HepLorentzVector pH,
double t1, double t2
) const{
ignore(p2out, p2in);
const double shat = p1in.invariantMass2(p2in);
assert(RHEJ::is_parton(id));
if(id != RHEJ::pid::gluon){
return 9./2.*shat*shat*C2qHqm(p1in,p1out,pH)/(t1*t2);
}
// gluon case
#ifdef RHEJ_BUILD_WITH_QCDLOOP
- if(!param.Hgg_setting().use_impact_factors){
+ if(!param_.Hgg_setting().use_impact_factors){
return C_A/C_F*1./(16*M_PI*M_PI)*t1/t2*MH2gq_outsideH(
p1out, p1in, p2out, p2in, pH,
- param.Hgg_setting().mt, param.Hgg_setting().include_bottom,
- param.Hgg_setting().mb
+ param_.Hgg_setting().mt, param_.Hgg_setting().include_bottom,
+ param_.Hgg_setting().mb
);
}
#endif
return 9./2.*shat*shat*(
C2gHgp(p1in,p1out,pH) + C2gHgm(p1in,p1out,pH)
)/(t1*t2);
}
double MatrixElement::tree_kin_Higgs_first(
std::array<Sparticle, 2> const & incoming,
std::vector<Sparticle> const & outgoing,
bool check_momenta
) const {
assert(outgoing.front().type == pid::Higgs);
const auto pH = to_HepLorentzVector(outgoing.front());
const auto partons = tag_extremal_jet_partons(
incoming,
std::vector<Sparticle>(begin(outgoing) + 1, end(outgoing)),
check_momenta
);
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto p1 = to_HepLorentzVector(partons.front());
const auto pn = to_HepLorentzVector(partons.back());
const auto q0 = pa - p1 - pH;
const double t1 = q0.m2();
const double t2 = (pn - pb).m2();
double wt = MH2_forwardH(
incoming[0].type, p1, pa, pn, pb, pH,
t1, t2
)*FKL_ladder_weight(
begin(partons) + 1, end(partons) - 1,
q0, pa, pb, p1, pn
);
for(auto const & inc: incoming){
if(inc.type != pid::gluon) wt *= C_F/C_A;
}
return wt;
}
double MatrixElement::tree_kin_Higgs_last(
std::array<Sparticle, 2> const & incoming,
std::vector<Sparticle> const & outgoing,
bool check_momenta
) const {
assert(outgoing.back().type == pid::Higgs);
const auto pH = to_HepLorentzVector(outgoing.back());
const auto partons = tag_extremal_jet_partons(
incoming,
std::vector<Sparticle>(begin(outgoing), end(outgoing) - 1),
check_momenta
);
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
auto p1 = to_HepLorentzVector(partons.front());
const auto pn = to_HepLorentzVector(partons.back());
auto q0 = pa - p1;
const double t1 = q0.m2();
const double t2 = (pn + pH - pb).m2();
double wt = MH2_forwardH(
incoming[1].type, pn, pb, p1, pa, pH,
t2, t1
)*FKL_ladder_weight(
begin(partons) + 1, end(partons) - 1,
q0, pa, pb, p1, pn
);
for(auto const & inc: incoming){
if(inc.type != pid::gluon) wt *= C_F/C_A;
}
return wt;
}
double MatrixElement::tree_kin_Higgs_between(
std::array<Sparticle, 2> const & incoming,
std::vector<Sparticle> const & outgoing,
bool check_momenta
) const {
const auto the_Higgs = std::find_if(
begin(outgoing), end(outgoing),
[](Sparticle const & s){ return s.type == pid::Higgs; }
);
assert(the_Higgs != end(outgoing));
const auto pH = to_HepLorentzVector(*the_Higgs);
std::vector<Sparticle> partons(begin(outgoing), the_Higgs);
partons.insert(end(partons), the_Higgs + 1, end(outgoing));
partons = tag_extremal_jet_partons(incoming, partons, check_momenta);
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
auto p1 = to_HepLorentzVector(
partons[has_unob_gluon(incoming, outgoing)?1:0]
);
auto pn = to_HepLorentzVector(
partons[partons.size() - (has_unof_gluon(incoming, outgoing)?2:1)]
);
auto first_after_Higgs = begin(partons) + (the_Higgs-begin(outgoing));
assert(
(first_after_Higgs == end(partons) && has_unob_gluon(incoming, outgoing))
|| first_after_Higgs->rapidity() >= the_Higgs->rapidity()
);
assert(
(first_after_Higgs == begin(partons) && has_unof_gluon(incoming, outgoing))
|| (first_after_Higgs-1)->rapidity() <= the_Higgs->rapidity()
);
// always treat the Higgs as if it were in between the extremal FKL partons
if(first_after_Higgs == begin(partons)) ++first_after_Higgs;
else if(first_after_Higgs == end(partons)) --first_after_Higgs;
// t-channel momentum before Higgs
auto qH = pa;
for(auto parton_it = begin(partons); parton_it != first_after_Higgs; ++parton_it){
qH -= to_HepLorentzVector(*parton_it);
}
auto q0 = pa - p1;
auto begin_ladder = begin(partons) + 1;
auto end_ladder = end(partons) - 1;
double current_factor;
if(has_unob_gluon(incoming, outgoing)){
current_factor = 9./2.*ME_Higgs_current_unob(
incoming[0].type, incoming[1].type,
pn, pb, to_HepLorentzVector(partons.front()), p1, pa, qH, qH - pH,
- param.Hgg_setting().mt, param.Hgg_setting().include_bottom, param.Hgg_setting().mb
+ param_.Hgg_setting().mt, param_.Hgg_setting().include_bottom, param_.Hgg_setting().mb
);
const auto p_unob = to_HepLorentzVector(partons.front());
q0 -= p_unob;
p1 += p_unob;
++begin_ladder;
}
else if(has_unof_gluon(incoming, outgoing)){
current_factor = 9./2.*ME_Higgs_current_unof(
incoming[0].type, incoming[1].type,
to_HepLorentzVector(partons.back()), pn, pb, p1, pa, qH, qH - pH,
- param.Hgg_setting().mt, param.Hgg_setting().include_bottom, param.Hgg_setting().mb
+ param_.Hgg_setting().mt, param_.Hgg_setting().include_bottom, param_.Hgg_setting().mb
);
pn += to_HepLorentzVector(partons.back());
--end_ladder;
}
else{
current_factor = ME_Higgs_current(
incoming[0].type, incoming[1].type,
pn, pb, p1, pa, qH, qH - pH,
- param.Hgg_setting().mt, param.Hgg_setting().include_bottom, param.Hgg_setting().mb
+ param_.Hgg_setting().mt, param_.Hgg_setting().include_bottom, param_.Hgg_setting().mb
);
}
const double ladder_factor = FKL_ladder_weight(
begin_ladder, first_after_Higgs,
q0, pa, pb, p1, pn
)*FKL_ladder_weight(
first_after_Higgs, end_ladder,
qH - pH, pa, pb, p1, pn
);
return current_factor*9./8.*ladder_factor;
}
double MatrixElement::tree_param_partons(
double alpha_s, double mur,
std::vector<Sparticle> const & partons
) const{
const double gs2 = 4.*M_PI*alpha_s;
double wt = std::pow(gs2, partons.size());
- if(param.log_corr()){
+ if(param_.log_corr()){
// use alpha_s(q_perp), evolved to mur
assert(partons.size() >= 2);
for(size_t i = 1; i < partons.size()-1; ++i){
wt *= 1 + alpha_s/(2*M_PI)*beta0*log(mur/partons[i].p.perp());
}
}
return wt;
}
double MatrixElement::tree_param(
double mur,
std::array<Sparticle, 2> const & incoming,
std::vector<Sparticle> const & outgoing
) const{
const double alpha_s = alpha_s_(mur);
if(has_unob_gluon(incoming, outgoing)){
return 4*M_PI*alpha_s*tree_param_partons(
alpha_s, mur, filter_partons({begin(outgoing) + 1, end(outgoing)})
);
}
if(has_unof_gluon(incoming, outgoing)){
return 4*M_PI*alpha_s*tree_param_partons(
alpha_s, mur, filter_partons({begin(outgoing), end(outgoing) - 1})
);
}
return tree_param_partons(alpha_s, mur, filter_partons(outgoing));
}
double MatrixElement::tree(
double mur,
std::array<Sparticle, 2> const & incoming,
std::vector<Sparticle> const & outgoing,
bool check_momenta
) const {
return tree_param(mur, incoming, outgoing)*tree_kin(
incoming, outgoing, check_momenta
);
}
} // namespace RHEJ
diff --git a/src/PhaseSpacePoint.cc b/src/PhaseSpacePoint.cc
index 1ecd7c3..c4cc22e 100644
--- a/src/PhaseSpacePoint.cc
+++ b/src/PhaseSpacePoint.cc
@@ -1,581 +1,581 @@
#include "RHEJ/PhaseSpacePoint.hh"
#include <random>
#include "RHEJ/Constants.hh"
#include "RHEJ/resummation_jet_momenta.hh"
#include "RHEJ/Jacobian.hh"
#include "RHEJ/RNGWrapper.hh"
#include "RHEJ/uno.hh"
#include "RHEJ/debug.hh"
#include "RHEJ/kinematics.hh"
namespace RHEJ{
namespace{
//generate Ranlux64Engine with fixed, predefined state
/*
* some (all?) of the Ranlux64Engine constructors leave fields
* uninitialised, invoking undefined behaviour. This can be
* circumvented by restoring the state from a file
*/
CLHEP::Ranlux64Engine gen_Ranlux64Engine(){
static const std::string state =
"9876\n"
"0.91280703978419097666\n"
"0.41606065829518357191\n"
"0.99156342622341142601\n"
"0.030922955274050423213\n"
"0.16206278421638486975\n"
"0.76151768001958330956\n"
"0.43765760066092695979\n"
"0.42904698253748563275\n"
"0.11476317525663759511\n"
"0.026620053590963976831\n"
"0.65953715764414511114\n"
"0.30136722624439826745\n"
"3.5527136788005009294e-15 4\n"
"1 202\n";
const std::string file = std::tmpnam(nullptr);
{
std::ofstream out{file};
out << state;
}
CLHEP::Ranlux64Engine result;
result.restoreStatus(file.c_str());
return result;
}
}
CLHEP::Ranlux64Engine PhaseSpacePoint::ran_{gen_Ranlux64Engine()};
void PhaseSpacePoint::reset_ranlux(std::string const & init_file){
reset_ranlux(init_file.c_str());
}
void PhaseSpacePoint::reset_ranlux(char const * init_file){
ran_.restoreStatus(init_file);
}
namespace {
constexpr int max_jet_user_idx = PhaseSpacePoint::ng_max;
bool is_nonjet_parton(fastjet::PseudoJet const & parton){
assert(parton.user_index() != -1);
return parton.user_index() > max_jet_user_idx;
}
bool is_jet_parton(fastjet::PseudoJet const & parton){
assert(parton.user_index() != -1);
return parton.user_index() <= max_jet_user_idx;
}
// user indices for partons with extremal rapidity
constexpr int unob_idx = -5;
constexpr int unof_idx = -4;
constexpr int backward_FKL_idx = -3;
constexpr int forward_FKL_idx = -2;
}
namespace{
double estimate_ng_mean(std::vector<fastjet::PseudoJet> const & Born_jets){
const double delta_y =
Born_jets.back().rapidity() - Born_jets.front().rapidity();
assert(delta_y > 0);
// Formula derived from fit in reversed HEJ intro paper
return 0.975052*delta_y;
}
}
std::vector<fastjet::PseudoJet> PhaseSpacePoint::cluster_jets(
std::vector<fastjet::PseudoJet> const & partons
) const{
- fastjet::ClusterSequence cs(partons, param.jet_def());
- return cs.inclusive_jets(param.jet_ptmin());
+ fastjet::ClusterSequence cs(partons, param_.jet_def());
+ return cs.inclusive_jets(param_.jet_ptmin());
}
bool PhaseSpacePoint::pass_resummation_cuts(
std::vector<fastjet::PseudoJet> const & jets
) const{
return cluster_jets(jets).size() == jets.size();
}
int PhaseSpacePoint::sample_ng(std::vector<fastjet::PseudoJet> const & Born_jets){
const double ng_mean = estimate_ng_mean(Born_jets);
std::poisson_distribution<int> dist(ng_mean);
RNGWrapper<CLHEP::Ranlux64Engine> rng{ran_};
const int ng = dist(rng);
assert(ng >= 0);
assert(ng < ng_max);
weight_ *= std::tgamma(ng + 1)*std::exp(ng_mean)*std::pow(ng_mean, -ng);
return ng;
}
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),
[](Sparticle 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);
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(outgoing_[new_idx].type == AWZH_boson->type);
decays_.emplace(new_idx, decay_it->second);
}
assert(std::is_sorted(begin(outgoing_), end(outgoing_), rapidity_less{}));
}
PhaseSpacePoint::PhaseSpacePoint(Event const & ev, Config const & conf):
unob_{has_unob_gluon(ev.incoming(), ev.outgoing())},
unof_{!unob_ && has_unof_gluon(ev.incoming(), ev.outgoing())},
- param(conf),
- splitter_{param.jet_R(), param.jet_def(), param.jet_ptmin(), ran_}
+ param_(conf),
+ splitter_{param_.jet_R(), param_.jet_def(), param_.jet_ptmin(), ran_}
{
weight_ = 1;
const auto Born_jets = sorted_by_rapidity(ev.jets());
const int ng = sample_ng(Born_jets);
weight_ /= std::tgamma(ng + 1);
const int ng_jets = sample_ng_jets(ng, Born_jets);
std::vector<fastjet::PseudoJet> out_partons = gen_non_jet(
- ng - ng_jets, CMINPT, param.jet_ptmin()
+ ng - ng_jets, CMINPT, param_.jet_ptmin()
);
{
const auto qperp = std::accumulate(
begin(out_partons), end(out_partons),
fastjet::PseudoJet{}
);
const auto jets = reshuffle(Born_jets, qperp);
if(weight_ == 0.) return;
if(! pass_resummation_cuts(jets)){
weight_ = 0.;
return;
}
std::vector<fastjet::PseudoJet> jet_partons = split(jets, ng_jets);
if(weight_ == 0.) return;
rescale_rapidities(
out_partons,
most_backward_FKL(jet_partons).rapidity(),
most_forward_FKL(jet_partons).rapidity()
);
if(! cluster_jets(out_partons).empty()){
weight_ = 0.;
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.;
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 & p: out_partons){
outgoing_.emplace_back(Sparticle{pid::gluon, std::move(p)});
}
most_backward_FKL(outgoing_).type = ev.incoming().front().type;
most_forward_FKL(outgoing_).type = ev.incoming().back().type;
copy_AWZH_boson_from(ev);
assert(!outgoing_.empty());
reconstruct_incoming(ev.incoming());
}
std::vector<fastjet::PseudoJet> PhaseSpacePoint::gen_non_jet(
int count, double ptmin, double ptmax
){
// heuristic parameters for pt sampling
const double ptpar = 1.3 + count/5.;
const double temp1 = atan((ptmax - ptmin)/ptpar);
std::vector<fastjet::PseudoJet> partons(count);
for(size_t i = 0; i < (size_t) count; ++i){
const double r1 = ran_.flat();
const double pt = ptmin + ptpar*tan(r1*temp1);
const double temp2 = cos(r1*temp1);
const double phi = 2*M_PI*ran_.flat();
weight_ *= 2.0*M_PI*pt*ptpar*temp1/(temp2*temp2);
// we don't know the allowed rapidity span yet,
// set a random value to be rescaled later on
const double y = 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 + RHEJ::PhaseSpacePoint::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 partons;
}
void PhaseSpacePoint::rescale_rapidities(
std::vector<fastjet::PseudoJet> & partons,
double ymin, double ymax
){
constexpr double ep = 1e-7;
for(auto & parton: partons){
assert(0 <= parton.rapidity() && parton.rapidity() <= 1);
const double dy = ymax - ymin - 2*ep;
const double y = ymin + ep + dy*parton.rapidity();
parton.reset_momentum_PtYPhiM(parton.pt(), y, parton.phi());
weight_ *= dy;
assert(ymin <= parton.rapidity() && parton.rapidity() <= ymax);
}
}
namespace {
template<typename T, typename... Rest>
auto min(T const & a, T const & b, Rest&&... r) {
using std::min;
return min(a, min(b, std::forward<Rest>(r)...));
}
}
double PhaseSpacePoint::probability_in_jet(
std::vector<fastjet::PseudoJet> 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_R();
+ const double R = param_.jet_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<fastjet::PseudoJet> const & Born_jets
){
RNGWrapper<CLHEP::Ranlux64Engine> rng{ran_};
const double p_J = probability_in_jet(Born_jets);
std::binomial_distribution<> bin_dist(ng, p_J);
const int ng_J = bin_dist(rng);
weight_ *= std::pow(p_J, -ng_J)*std::pow(1 - p_J, ng_J - ng);
return ng_J;
}
std::vector<fastjet::PseudoJet>
PhaseSpacePoint::reshuffle(
std::vector<fastjet::PseudoJet> const & Born_jets,
fastjet::PseudoJet const & q
){
if(q == fastjet::PseudoJet{0, 0, 0, 0}) return Born_jets;
std::vector<fastjet::PseudoJet> jets = resummation_jet_momenta(Born_jets, q);
if(jets.empty()){
weight_ = 0;
return {};
}
// transform delta functions to integration over resummation momenta
weight_ /= Jacobian(jets, q);
return jets;
}
std::vector<int> PhaseSpacePoint::distribute_jet_partons(
int ng_jets, std::vector<fastjet::PseudoJet> const & jets
){
size_t first_valid_jet = 0;
size_t num_valid_jets = jets.size();
- const double R_eff = 5./3.*param.jet_R();
+ const double R_eff = 5./3.*param_.jet_R();
// if there is an unordered jet too far away from the FKL jets
// then extra gluon constituents of the unordered jet would
// violate the FKL rapidity ordering
if(unob_ && jets[0].delta_R(jets[1]) > R_eff){
++first_valid_jet;
--num_valid_jets;
}
else if(unof_ && jets[jets.size()-1].delta_R(jets[jets.size()-2]) > R_eff){
--num_valid_jets;
}
std::vector<int> np(jets.size(), 1);
for(int i = 0; i < ng_jets; ++i){
++np[first_valid_jet + ran_.flat() * num_valid_jets];
}
weight_ *= std::pow(num_valid_jets, ng_jets);
return np;
}
#ifndef NDEBUG
namespace{
bool tagged_FKL_backward(
std::vector<fastjet::PseudoJet> const & jet_partons
){
return std::find_if(
begin(jet_partons), end(jet_partons),
[](fastjet::PseudoJet const & p){
return p.user_index() == backward_FKL_idx;
}
) != end(jet_partons);
}
bool tagged_FKL_forward(
std::vector<fastjet::PseudoJet> 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() == forward_FKL_idx;
}
) != jet_partons.rend();
}
bool tagged_FKL_extremal(
std::vector<fastjet::PseudoJet> const & jet_partons
){
return tagged_FKL_backward(jet_partons) && tagged_FKL_forward(jet_partons);
}
} // namespace anonymous
#endif
std::vector<fastjet::PseudoJet> PhaseSpacePoint::split(
std::vector<fastjet::PseudoJet> const & jets,
int ng_jets
){
return split(jets, distribute_jet_partons(ng_jets, jets));
}
bool PhaseSpacePoint::pass_extremal_cuts(
fastjet::PseudoJet const & ext_parton,
fastjet::PseudoJet const & jet
) const{
- if(ext_parton.pt() < param.min_ext_pt()) return false;
- return (ext_parton - jet).pt()/jet.pt() < param.max_ext_soft_pt_fraction();
+ if(ext_parton.pt() < param_.min_ext_pt()) return false;
+ return (ext_parton - jet).pt()/jet.pt() < param_.max_ext_soft_pt_fraction();
}
std::vector<fastjet::PseudoJet> PhaseSpacePoint::split(
std::vector<fastjet::PseudoJet> const & jets,
std::vector<int> const & np
){
assert(! jets.empty());
assert(jets.size() == np.size());
assert(pass_resummation_cuts(jets));
const size_t most_backward_FKL_idx = 0 + unob_;
const size_t most_forward_FKL_idx = jets.size() - 1 - unof_;
std::vector<fastjet::PseudoJet> jet_partons;
// randomly distribute jet gluons among jets
for(size_t i = 0; i < jets.size(); ++i){
weight_ *= splitter_.Split(jets[i], np[i]);
if(weight_ == 0) return {};
assert(
std::all_of(
begin(splitter_.get_jcons()), end(splitter_.get_jcons()),
is_jet_parton
)
);
const auto first_new_parton = jet_partons.insert(
end(jet_partons),
begin(splitter_.get_jcons()), end(splitter_.get_jcons())
);
// mark uno and extremal FKL emissions here so we can check
// their position once all emissions are generated
auto extremal = end(jet_partons);
if((unob_ && i == 0) || i == most_backward_FKL_idx){
// unordered or FKL backward emission
extremal = std::min_element(
first_new_parton, end(jet_partons), rapidity_less{}
);
extremal->set_user_index(
(i == most_backward_FKL_idx)?backward_FKL_idx:unob_idx
);
}
else if((unof_ && i == jets.size() - 1) || i == most_forward_FKL_idx){
// unordered or FKL forward emission
extremal = std::max_element(
first_new_parton, end(jet_partons), rapidity_less{}
);
extremal->set_user_index(
(i == most_forward_FKL_idx)?forward_FKL_idx:unof_idx
);
}
if(
extremal != end(jet_partons)
&& !pass_extremal_cuts(*extremal, jets[i])
){
weight_ = 0;
return {};
}
}
assert(tagged_FKL_extremal(jet_partons));
std::sort(begin(jet_partons), end(jet_partons), rapidity_less{});
if(
!extremal_ok(jet_partons)
|| !split_preserved_jets(jets, jet_partons)
){
weight_ = 0.;
return {};
}
return jet_partons;
}
bool PhaseSpacePoint::extremal_ok(
std::vector<fastjet::PseudoJet> const & partons
) const{
assert(std::is_sorted(begin(partons), end(partons), rapidity_less{}));
if(unob_ && partons.front().user_index() != unob_idx) return false;
if(unof_ && partons.back().user_index() != unof_idx) return false;
return
most_backward_FKL(partons).user_index() == backward_FKL_idx
&& most_forward_FKL(partons).user_index() == forward_FKL_idx;
}
bool PhaseSpacePoint::split_preserved_jets(
std::vector<fastjet::PseudoJet> const & jets,
std::vector<fastjet::PseudoJet> const & jet_partons
) const{
assert(std::is_sorted(begin(jets), end(jets), rapidity_less{}));
const auto split_jets = sorted_by_rapidity(cluster_jets(jet_partons));
// this can happen if two overlapping jets
// are both split into more than one parton
if(split_jets.size() != jets.size()) return false;
for(size_t i = 0; i < split_jets.size(); ++i){
// this can happen if there are two overlapping jets
// and a parton is assigned to the "wrong" jet
if(!nearby_ep(jets[i].rapidity(), split_jets[i].rapidity(), 1e-2)){
return false;
}
}
return true;
}
template<class Particle>
Particle const & PhaseSpacePoint::most_backward_FKL(
std::vector<Particle> const & partons
) const{
return partons[0 + unob_];
}
template<class Particle>
Particle const & PhaseSpacePoint::most_forward_FKL(
std::vector<Particle> const & partons
) const{
const size_t idx = partons.size() - 1 - unof_;
assert(idx < partons.size());
return partons[idx];
}
template<class Particle>
Particle & PhaseSpacePoint::most_backward_FKL(
std::vector<Particle> & partons
) const{
return partons[0 + unob_];
}
template<class Particle>
Particle & PhaseSpacePoint::most_forward_FKL(
std::vector<Particle> & partons
) const{
const size_t idx = partons.size() - 1 - unof_;
assert(idx < partons.size());
return partons[idx];
}
namespace{
bool contains_idx(
fastjet::PseudoJet const & jet, fastjet::PseudoJet const & parton
){
auto const & constituents = jet.constituents();
const int idx = parton.user_index();
return std::find_if(
begin(constituents), end(constituents),
[idx](fastjet::PseudoJet const & con){return con.user_index() == idx;}
) != end(constituents);
}
}
/**
* final jet test:
* - number of jets must match Born kinematics
* - no partons designated as nonjet may end up inside jets
* - all other outgoing partons *must* end up inside jets
* - the extremal (in rapidity) partons must be inside the extremal jets
* - rapidities must be the same (by construction)
*/
bool PhaseSpacePoint::jets_ok(
std::vector<fastjet::PseudoJet> const & Born_jets,
std::vector<fastjet::PseudoJet> const & partons
) const{
- fastjet::ClusterSequence cs(partons, param.jet_def());
- const auto jets = sorted_by_rapidity(cs.inclusive_jets(param.jet_ptmin()));
+ fastjet::ClusterSequence cs(partons, param_.jet_def());
+ const auto jets = sorted_by_rapidity(cs.inclusive_jets(param_.jet_ptmin()));
if(jets.size() != Born_jets.size()) return false;
int in_jet = 0;
for(size_t i = 0; i < jets.size(); ++i){
assert(jets[i].has_constituents());
for(auto && parton: jets[i].constituents()){
if(is_nonjet_parton(parton)) return false;
}
in_jet += jets[i].constituents().size();
}
const int expect_in_jet = std::count_if(
partons.cbegin(), partons.cend(), is_jet_parton
);
if(in_jet != expect_in_jet) return false;
// note that PseudoJet::contains does not work here
if(! (
contains_idx(most_backward_FKL(jets), most_backward_FKL(partons))
&& contains_idx(most_forward_FKL(jets), most_forward_FKL(partons))
)) return false;
if(unob_ && !contains_idx(jets.front(), partons.front())) return false;
if(unof_ && !contains_idx(jets.back(), partons.back())) return false;
for(size_t i = 0; i < jets.size(); ++i){
assert(nearby_ep(jets[i].rapidity(), Born_jets[i].rapidity(), 1e-2));
}
return true;
}
namespace{
}
void PhaseSpacePoint::reconstruct_incoming(
std::array<Sparticle, 2> const & Born_incoming
){
std::tie(incoming_[0].p, incoming_[1].p) = incoming_momenta(outgoing_);
for(size_t i = 0; i < incoming_.size(); ++i){
incoming_[i].type = Born_incoming[i].type;
}
assert(momentum_conserved(1e-7));
}
double PhaseSpacePoint::phase_space_normalisation(
int num_Born_jets, int num_out_partons
) const{
return pow(16*pow(M_PI,3), num_Born_jets - num_out_partons);
}
bool PhaseSpacePoint::momentum_conserved(double ep) const{
fastjet::PseudoJet diff;
for(auto const & in: incoming()) diff += in.p;
for(auto const & out: outgoing()) diff -= out.p;
return nearby_ep(diff, fastjet::PseudoJet{}, ep);
}
}
diff --git a/src/main.cc b/src/main.cc
index 050924c..90fa1c6 100644
--- a/src/main.cc
+++ b/src/main.cc
@@ -1,162 +1,159 @@
/**
* Name: main.cc
* Authors: Tuomas Hapola, Andreas Maier <andreas.maier@durham.ac.uk>
*
*/
#include <fstream>
#include <algorithm>
#include <memory>
#include <chrono>
#include <iostream>
#include "yaml-cpp/yaml.h"
#include "RHEJ/CombinedEventWriter.hh"
#include "RHEJ/config.hh"
#include "RHEJ/EventReweighter.hh"
#include "RHEJ/get_analysis.hh"
#include "LHEF/LHEF.h"
#include "RHEJ/utility.hh"
#include "RHEJ/Version.hh"
#include "RHEJ/stream.hh"
#include "RHEJ/YAMLreader.hh"
int event_number(std::string const & record){
size_t start = record.rfind("Number of Events");
start = record.find_first_of("123456789", start);
if(start == std::string::npos) {
throw std::invalid_argument("no event number record found");
}
const size_t end = record.find_first_not_of("0123456789", start);
return std::stoi(record.substr(start, end - start));
}
RHEJ::Config load_config(char const * filename){
try{
return RHEJ::load_config(filename);
}
catch(std::exception const & exc){
std::cerr << "Error: " << exc.what() << '\n';
std::exit(EXIT_FAILURE);
}
}
std::unique_ptr<RHEJ::Analysis> get_analysis(
YAML::Node const & parameters
){
try{
return RHEJ::get_analysis(parameters);
}
catch(std::exception const & exc){
std::cerr << "Failed to load analysis: " << exc.what() << '\n';
std::exit(EXIT_FAILURE);
}
}
-std::string print_time(const time_t time){
+std::string time_to_string(const time_t time){
char s[1000];
struct tm * p = localtime(&time);
strftime(s, 1000, "%a %b %d %Y %H:%M:%S", p);
return s;
}
int main(int argn, char** argv) {
using clock = std::chrono::system_clock;
if (argn < 3) {
std::cerr << "\n# Usage:\n."<< argv[0] <<" config_file input_file\n\n";
return EXIT_FAILURE;
}
const auto start_time = clock::now();
{
- // char tmp_s[30];
- // const auto tmp_time = clock::to_time_t(start_time);
- // strftime(tmp_s, 30, "%a %b %d %Y %H:%M:%S", localtime(&(tmp_time)));
std::cout << "Starting " << RHEJ::package_name_full() << " ("
- << print_time(clock::to_time_t(start_time)) << ")" << std::endl;
+ << time_to_string(clock::to_time_t(start_time)) << ")" << std::endl;
}
// read configuration
const RHEJ::Config config = load_config(argv[1]);
RHEJ::istream in{argv[2]};
std::unique_ptr<RHEJ::Analysis> analysis = get_analysis(
config.analysis_parameters
);
assert(analysis != nullptr);
LHEF::Reader reader{in};
RHEJ::CombinedEventWriter writer{config, reader.heprup};
double global_reweight = 1.;
int max_events = std::numeric_limits<int>::max();
if(argn > 3){
max_events = std::stoi(argv[3]);
const int input_events = event_number(reader.headerBlock);
global_reweight = input_events/static_cast<double>(max_events);
std::cout << "Processing " << max_events
<< " out of " << input_events << " events\n";
}
RHEJ::EventReweighter rhej{
reader.heprup,
config
};
if(config.RanLux_init){
RHEJ::PhaseSpacePoint::reset_ranlux(*config.RanLux_init);
}
int nevent = 0;
std::array<int, RHEJ::event_type::last_type + 1>
nevent_type{0}, nfailed_type{0};
// Loop over the events in the inputfile
while(reader.readEvent()){
// reweight events so that the total cross section is conserved
reader.hepeup.setWeight(0, global_reweight * reader.hepeup.weight());
if(nevent == max_events) break;
++nevent;
if (nevent % 10000 == 0){
std::cout << "Passed " << nevent << " events ("
- << print_time(clock::to_time_t(clock::now())) << ")"<< std::endl;
+ << time_to_string(clock::to_time_t(clock::now())) << ")"<< std::endl;
}
// calculate rHEJ weight
RHEJ::Event FO_event{
RHEJ::UnclusteredEvent{reader.hepeup},
config.fixed_order_jets.def, config.fixed_order_jets.min_pt,
};
auto resummed_events = rhej.reweight(
FO_event,
config.trials, config.scale_gen
);
++nevent_type[FO_event.type()];
if(resummed_events.empty()) ++nfailed_type[FO_event.type()];
for(auto const & ev: resummed_events){
//TODO: move pass_cuts to after phase space point generation
if(analysis->pass_cuts(ev, FO_event)){
analysis->fill(ev, FO_event);
writer.write(ev);
}
}
} // main event loop
using namespace RHEJ::event_type;
std::cout<< "Events processed: " << nevent << '\n';
for(size_t ev_type = first_type; ev_type <= last_type; ++ev_type){
std::cout << '\t' << names[ev_type] << ": " << nevent_type[ev_type]
<< ", failed to reconstruct " << nfailed_type[ev_type]
<< '\n';
}
std::chrono::duration<double> run_time = (clock::now() - start_time);
std::cout << "Finished " << RHEJ::package_name() << " at "
- << print_time(clock::to_time_t(clock::now()))
+ << time_to_string(clock::to_time_t(clock::now()))
<< "\n=> Runtime: " << run_time.count() << " sec ("
<< nevent/run_time.count() << " Events/sec).\n";
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Jan 21, 1:47 AM (1 d, 11 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4243532
Default Alt Text
(106 KB)

Event Timeline