Value with uncertainty object

CepGen/Core/Generator.cpp

@@ -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();
   }
 
@@ -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() {
@@ -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); }
@@ -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();

CepGen/EventFilter/EventExporter.h

@@ -25,6 +25,7 @@
 
 namespace cepgen {
   class Event;
+  class Value;
   /**
    * \brief Output format handler for events export
    * \author Laurent Forthomme <laurent.forthomme@cern.ch>
@@ -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; }
 
@@ -51,4 +52,3 @@ namespace cepgen {
 }  // namespace cepgen
 
 #endif
-

CepGen/EventFilter/EventHarvester.h

@@ -20,6 +20,7 @@
 
 #include "CepGen/EventFilter/EventExporter.h"
 #include "CepGen/Utils/Histogram.h"
+#include "CepGen/Utils/Value.h"
 
 namespace cepgen {
   class Event;
@@ -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:
@@ -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
@@ -76,4 +77,3 @@ namespace cepgen {
     std::vector<Hist2DInfo> hists2d_;
   };
 }  // namespace cepgen
-

CepGen/EventFilter/EventModifier.h

@@ -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
@@ -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

CepGen/Generator.h

@@ -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
@@ -19,6 +19,7 @@
 #include <memory>
 
 #include "CepGen/Event/Event.h"
+#include "CepGen/Utils/Value.h"
 
 ////////////////////////////////////////////////////////////////////////////////
 
@@ -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.
@@ -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
 

CepGen/Integration/Integrator.cpp

@@ -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);

CepGen/Integration/Integrator.h

@@ -23,6 +23,7 @@
 #include <vector>
 
 #include "CepGen/Modules/NamedModule.h"
+#include "CepGen/Utils/Value.h"
 
 namespace cepgen {
   class Integrand;
@@ -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_;
   };

CepGen/Integration/IntegratorMISER.cpp

@@ -30,7 +30,7 @@ namespace cepgen {
 
     static ParametersDescription description();
 
-    void integrate(Integrand&, double&, double&) override;
+    Value integrate(Integrand&) override;
 
   private:
     int ncvg_;
@@ -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_;
@@ -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],
@@ -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() {

CepGen/Integration/IntegratorPlain.cpp

@@ -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_;
@@ -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],
@@ -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() {

CepGen/Integration/IntegratorVegas.cpp

@@ -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 };
 
@@ -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
@@ -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],
@@ -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) {