Skip to content

Commit

Permalink
Value with uncertainty object (#29)
Browse files Browse the repository at this point in the history 
* Core: new Value object introduced

* Integrator: using Value wherever possible

* Utils: using Value in histograms/graphs

* Tests: using Value object in processes cross sections comparator

* Utils: templated arithmetic operators for values

* Core: removed a batch of unnecessary explicit casts of value object to double
  • Loading branch information
@forthommel
forthommel committed Sep 13, 2023
1 parent 111c71a commit 031e355
Show file tree
Hide file tree
Showing 53 changed files with 367 additions and 266 deletions.
41 changes: 23 additions & 18 deletions CepGen/Core/Generator.cpp
View file
Expand Up @@ -64,7 +64,7 @@ namespace cepgen {

worker_->setRuntimeParameters(const_cast<const Parameters*>(parameters_.get()));
worker_->setIntegrator(integrator_.get());
result_ = result_error_ = -1.;
xsect_ = Value{-1., -1.};
parameters_->prepareRun();
}

Expand All @@ -86,23 +86,28 @@ namespace cepgen {
}

void Generator::computeXsection(double& cross_section, double& err) {
const auto xsec = computeXsection();
cross_section = xsec;
err = xsec.uncertainty();
}

Value Generator::computeXsection() {
CG_INFO("Generator") << "Starting the computation of the process cross-section.";

integrate(); // run is cleared here

cross_section = result_;
err = result_error_;

if (cross_section < 1.e-2)
CG_INFO("Generator") << "Total cross section: " << cross_section * 1.e3 << " +/- " << err * 1.e3 << " fb.";
else if (cross_section < 0.5e3)
CG_INFO("Generator") << "Total cross section: " << cross_section << " +/- " << err << " pb.";
else if (cross_section < 0.5e6)
CG_INFO("Generator") << "Total cross section: " << cross_section * 1.e-3 << " +/- " << err * 1.e-3 << " nb.";
else if (cross_section < 0.5e9)
CG_INFO("Generator") << "Total cross section: " << cross_section * 1.e-6 << " +/- " << err * 1.e-6 << " µb.";
if (xsect_ < 1.e-2)
CG_INFO("Generator") << "Total cross section: " << xsect_ * 1.e3 << " fb.";
else if (xsect_ < 0.5e3)
CG_INFO("Generator") << "Total cross section: " << xsect_ << " pb.";
else if (xsect_ < 0.5e6)
CG_INFO("Generator") << "Total cross section: " << xsect_ * 1.e-3 << " nb.";
else if (xsect_ < 0.5e9)
CG_INFO("Generator") << "Total cross section: " << xsect_ * 1.e-6 << " µb.";
else
CG_INFO("Generator") << "Total cross section: " << cross_section * 1.e-9 << " +/- " << err * 1.e-9 << " mb.";
CG_INFO("Generator") << "Total cross section: " << xsect_ * 1.e-9 << " mb.";

return xsect_;
}

void Generator::resetIntegrator() {
Expand Down Expand Up @@ -137,16 +142,16 @@ namespace cepgen {
if (!integrator_)
throw CG_FATAL("Generator:integrate") << "No integrator object was declared for the generator!";

integrator_->integrate(worker_->integrand(), result_, result_error_);
xsect_ = integrator_->integrate(worker_->integrand());

CG_DEBUG("Generator:integrate") << "Computed cross section: (" << result_ << " +- " << result_error_ << ") pb.";
CG_DEBUG("Generator:integrate") << "Computed cross section: (" << xsect_ << ") pb.";

// now that the cross section has been computed, feed it to the event modification algorithms...
for (auto& mod : parameters_->eventModifiersSequence())
mod->setCrossSection(result_, result_error_);
mod->setCrossSection(xsect_);
// ...and to the event storage algorithms
for (auto& mod : parameters_->eventExportersSequence())
mod->setCrossSection(result_, result_error_);
mod->setCrossSection(xsect_);
}

const Event& Generator::generateOneEvent(Event::callback callback) { return next(callback); }
Expand Down Expand Up @@ -188,7 +193,7 @@ namespace cepgen {
//--- if invalid argument, retrieve from runtime parameters
if (num_events < 1) {
if (parameters_->generation().targetLuminosity() > 0.) {
num_events = std::ceil(parameters_->generation().targetLuminosity() * result_);
num_events = std::ceil(parameters_->generation().targetLuminosity() * xsect_);
CG_INFO("Generator") << "Target luminosity: " << parameters_->generation().targetLuminosity() << " pb-1.";
} else
num_events = parameters_->generation().maxGen();
Expand Down
6 changes: 3 additions & 3 deletions CepGen/EventFilter/EventExporter.h
View file
Expand Up @@ -25,6 +25,7 @@

namespace cepgen {
class Event;
class Value;
/**
* \brief Output format handler for events export
* \author Laurent Forthomme <laurent.forthomme@cern.ch>
Expand All @@ -34,8 +35,8 @@ namespace cepgen {
public:
explicit EventExporter(const ParametersList&);

/// Set the process cross section and its associated error
virtual void setCrossSection(double /*cross_section*/, double /*err_cross_section*/) {}
/// Specify the process cross section and uncertainty, in pb
virtual void setCrossSection(const Value&) {}
/// Set the event number
void setEventNumber(unsigned long long ev_id) { event_num_ = ev_id; }

Expand All @@ -51,4 +52,3 @@ namespace cepgen {
} // namespace cepgen

#endif

6 changes: 3 additions & 3 deletions CepGen/EventFilter/EventHarvester.h
View file
Expand Up @@ -20,6 +20,7 @@

#include "CepGen/EventFilter/EventExporter.h"
#include "CepGen/Utils/Histogram.h"
#include "CepGen/Utils/Value.h"

namespace cepgen {
class Event;
Expand All @@ -40,7 +41,7 @@ namespace cepgen {
static ParametersDescription description();

void initialise() override;
void setCrossSection(double cross_section, double) override { cross_section_ = cross_section; }
void setCrossSection(const Value& cross_section) override { cross_section_ = cross_section; }
void operator<<(const Event&) override;

private:
Expand All @@ -52,7 +53,7 @@ namespace cepgen {
std::ofstream file_;
std::unique_ptr<utils::Drawer> drawer_;

double cross_section_{1.};
Value cross_section_{1., 0.};
unsigned long num_evts_{0ul};

/// Name of the physics process
Expand All @@ -76,4 +77,3 @@ namespace cepgen {
std::vector<Hist2DInfo> hists2d_;
};
} // namespace cepgen

3 changes: 2 additions & 1 deletion CepGen/EventFilter/EventModifier.h
View file
Expand Up @@ -25,6 +25,7 @@
#include "CepGen/EventFilter/EventHandler.h"

namespace cepgen {
class Value;
/// Class template to interface (external/internal) events modification algorithms
/// \author Laurent Forthomme <laurent.forthomme@cern.ch>
/// \date July 2019
Expand Down Expand Up @@ -55,7 +56,7 @@ namespace cepgen {
*/
virtual bool run(Event& ev, double& weight, bool full = true) = 0;
/// Specify the process cross section and uncertainty, in pb
virtual void setCrossSection(double, double) {}
virtual void setCrossSection(const Value&) {}

protected:
/// Random numbers generator seed fed to the algorithm
Expand Down
26 changes: 12 additions & 14 deletions CepGen/Generator.h
View file
@@ -1,6 +1,6 @@
/*
* CepGen: a central exclusive processes event generator
* Copyright (C) 2013-2022 Laurent Forthomme
* Copyright (C) 2013-2023 Laurent Forthomme
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
Expand All @@ -19,6 +19,7 @@
#include <memory>

#include "CepGen/Event/Event.h"
#include "CepGen/Utils/Value.h"

////////////////////////////////////////////////////////////////////////////////

Expand Down Expand Up @@ -112,18 +113,17 @@ namespace cepgen {
void clearRun();
/// Integrate the functional over the whole phase space
void integrate();
/**
* Compute the cross section for the run parameters defined by this object.
* This returns the cross section as well as the absolute error computed along.
* \brief Compute the cross-section for the given process
* \param[out] xsec The computed cross-section, in pb
* \param[out] err The absolute integration error on the computed cross-section, in pb
*/
void computeXsection(double& cross_section, double& err);
/// Compute the cross section for the run parameters
/// \return The computed cross-section and uncertainty, in pb
Value computeXsection();
/// Compute the cross section for the run parameters
/// \param[out] xsec The computed cross-section, in pb
/// \param[out] err The absolute integration error on the computed cross-section, in pb
[[deprecated("Please use the parameters-less version")]] void computeXsection(double& cross_section, double& err);
/// Last cross section computed by the generator
double crossSection() const { return result_; }
double crossSection() const { return xsect_; }
/// Last error on the cross section computed by the generator
double crossSectionError() const { return result_error_; }
double crossSectionError() const { return xsect_.uncertainty(); }

// void terminate();
/// \deprecated Replaced by the generic method \a generate.
Expand Down Expand Up @@ -152,9 +152,7 @@ namespace cepgen {
/// Has the event generator already been initialised?
bool initialised_{false};
/// Cross section value computed at the last integration
double result_{-1.};
/// Error on the cross section as computed in the last integration
double result_error_{-1.};
Value xsect_{-1., -1.};
};
} // namespace cepgen

Expand Down
18 changes: 8 additions & 10 deletions CepGen/Integration/Integrator.cpp
View file
Expand Up @@ -53,23 +53,21 @@ namespace cepgen {

double Integrator::uniform(double min, double max) const { return min + (max - min) * rnd_(rnd_gen_); }

double Integrator::integrate(Integrand& integrand) {
Value Integrator::integrate(Integrand& integrand) {
if (limits_.size() != integrand.size())
limits_ = std::vector<Limits>(integrand.size(), Limits{0., 1.});
double result, tmp;
integrate(integrand, result, tmp);
return result;
return integrate(integrand);
}

double Integrator::integrate(const std::function<double(const std::vector<double>&)>& func,
const ParametersList& params,
size_t num_vars) {
Value Integrator::integrate(const std::function<double(const std::vector<double>&)>& func,
const ParametersList& params,
size_t num_vars) {
return integrate(func, params, std::vector<Limits>(num_vars, Limits{0., 1.}));
}

double Integrator::integrate(const std::function<double(const std::vector<double>&)>& func,
const ParametersList& params,
const std::vector<Limits>& limits) {
Value Integrator::integrate(const std::function<double(const std::vector<double>&)>& func,
const ParametersList& params,
const std::vector<Limits>& limits) {
auto integr = IntegratorFactory::get().build(params);
integr->setLimits(limits);
auto integrand = FunctionIntegrand(limits.size(), func);
Expand Down
17 changes: 6 additions & 11 deletions CepGen/Integration/Integrator.h
View file
Expand Up @@ -23,6 +23,7 @@
#include <vector>

#include "CepGen/Modules/NamedModule.h"
#include "CepGen/Utils/Value.h"

namespace cepgen {
class Integrand;
Expand All @@ -46,24 +47,18 @@ namespace cepgen {

/// Perform the multidimensional Monte Carlo integration
/// \param[out] result integral computed over the full phase space
/// \param[out] abserr uncertainty associated to the computed integral
virtual void integrate(Integrand&, double& result, double& abserr) = 0;
/// Perform an integration with no use of the numerical error
/// \return the integral computed over the full phase space
double integrate(Integrand&);
virtual Value integrate(Integrand& result) = 0;
/// Perform an integration with a given functional and a given set of parameters
static double integrate(const std::function<double(const std::vector<double>&)>&, const ParametersList&, size_t);
static Value integrate(const std::function<double(const std::vector<double>&)>&, const ParametersList&, size_t);
/// Perform an integration with a given functional and a given set of parameters
static double integrate(const std::function<double(const std::vector<double>&)>&,
const ParametersList&,
const std::vector<Limits>&);
static Value integrate(const std::function<double(const std::vector<double>&)>&,
const ParametersList&,
const std::vector<Limits>&);

protected:
const unsigned long seed_; ///< Random number generator seed
int verbosity_; ///< Integrator verbosity
std::vector<Limits> limits_; ///< List of per-variable integration limits
double result_{0.}; ///< Result of the last integration
double err_result_{0.}; ///< Standard deviation for the last integration
mutable std::default_random_engine rnd_gen_;
mutable std::uniform_real_distribution<double> rnd_;
};
Expand Down
8 changes: 4 additions & 4 deletions CepGen/Integration/IntegratorMISER.cpp
View file
Expand Up @@ -30,7 +30,7 @@ namespace cepgen {

static ParametersDescription description();

void integrate(Integrand&, double&, double&) override;
Value integrate(Integrand&) override;

private:
int ncvg_;
Expand All @@ -45,7 +45,7 @@ namespace cepgen {
IntegratorMISER::IntegratorMISER(const ParametersList& params)
: IntegratorGSL(params), ncvg_(steer<int>("numFunctionCalls")) {}

void IntegratorMISER::integrate(Integrand& integrand, double& result, double& abserr) {
Value IntegratorMISER::integrate(Integrand& integrand) {
setIntegrand(integrand);
miser_state_.reset(gsl_monte_miser_alloc(function_->dim));
miser_state_->verbose = verbosity_;
Expand All @@ -65,6 +65,7 @@ namespace cepgen {
<< "Dither: " << miser_params_.dither << ".";

// launch the full integration
double result, abserr;
int res = gsl_monte_miser_integrate(function_.get(),
&xlow_[0],
&xhigh_[0],
Expand All @@ -79,8 +80,7 @@ namespace cepgen {
throw CG_FATAL("Integrator:integrate") << "Error while performing the integration!\n\t"
<< "GSL error: " << gsl_strerror(res) << ".";

result_ = result;
err_result_ = abserr;
return Value{result, abserr};
}

ParametersDescription IntegratorMISER::description() {
Expand Down
8 changes: 4 additions & 4 deletions CepGen/Integration/IntegratorPlain.cpp
View file
Expand Up @@ -30,7 +30,7 @@ namespace cepgen {

static ParametersDescription description();

void integrate(Integrand&, double& result, double& abserr) override;
Value integrate(Integrand&) override;

private:
const int ncvg_;
Expand All @@ -39,12 +39,13 @@ namespace cepgen {
IntegratorPlain::IntegratorPlain(const ParametersList& params)
: IntegratorGSL(params), ncvg_(steer<int>("numFunctionCalls")) {}

void IntegratorPlain::integrate(Integrand& integrand, double& result, double& abserr) {
Value IntegratorPlain::integrate(Integrand& integrand) {
setIntegrand(integrand);

//--- launch integration
std::unique_ptr<gsl_monte_plain_state, void (*)(gsl_monte_plain_state*)> pln_state(
gsl_monte_plain_alloc(function_->dim), gsl_monte_plain_free);
double result, abserr;
int res = gsl_monte_plain_integrate(function_.get(),
&xlow_[0],
&xhigh_[0],
Expand All @@ -59,8 +60,7 @@ namespace cepgen {
throw CG_FATAL("Integrator:integrate") << "Error while performing the integration!\n\t"
<< "GSL error: " << gsl_strerror(res) << ".";

result_ = result;
err_result_ = abserr;
return Value{result, abserr};
}

ParametersDescription IntegratorPlain::description() {
Expand Down
8 changes: 4 additions & 4 deletions CepGen/Integration/IntegratorVegas.cpp
View file
Expand Up @@ -32,7 +32,7 @@ namespace cepgen {

static ParametersDescription description();

void integrate(Integrand&, double&, double&) override;
Value integrate(Integrand&) override;

enum class Mode { importance = 1, importanceOnly = 0, stratified = -1 };

Expand Down Expand Up @@ -68,7 +68,7 @@ namespace cepgen {
verbosity_ = steer<int>("verbose"); // supersede the parent default verbosity level
}

void IntegratorVegas::integrate(Integrand& integrand, double& result, double& abserr) {
Value IntegratorVegas::integrate(Integrand& integrand) {
setIntegrand(integrand);

//--- start by preparing the grid/state
Expand Down Expand Up @@ -106,6 +106,7 @@ namespace cepgen {

// integration phase
unsigned short it_chisq = 0;
double result, abserr;
do {
int res = gsl_monte_vegas_integrate(function_.get(),
&xlow_[0],
Expand Down Expand Up @@ -133,8 +134,7 @@ namespace cepgen {
<< "and generated " << vegas_state_->bins_max << " bins.\n\t"
<< "Integration volume: " << vegas_state_->vol << ".";

result_ = result;
err_result_ = abserr;
return Value{result, abserr};
}

void IntegratorVegas::warmup(size_t ncall) {
Expand Down

0 comments on commit 031e355

Please sign in to comment.