From 268a5499f010e2897c30e0cc58f7a1ab71080541 Mon Sep 17 00:00:00 2001 From: Laurent Forthomme Date: Sat, 10 Feb 2024 17:25:59 +0100 Subject: [PATCH 1/7] Addition of a ROOT event generator instance --- .../event_generator/EventGeneratorROOT.h | 102 +++++++++ .../event_generator/EventGeneratorROOT.cpp | 210 ++++++++++++++++++ 2 files changed, 312 insertions(+) create mode 100644 include/modules/event_generator/EventGeneratorROOT.h create mode 100644 src/modules/event_generator/EventGeneratorROOT.cpp diff --git a/include/modules/event_generator/EventGeneratorROOT.h b/include/modules/event_generator/EventGeneratorROOT.h new file mode 100644 index 0000000..320c718 --- /dev/null +++ b/include/modules/event_generator/EventGeneratorROOT.h @@ -0,0 +1,102 @@ +/* + * EventGeneratorROOT.h + * + * Created on: Feb 2, 2024 + * Author: Laurent Forthomme (AGH) + */ + +#ifndef MODULES_EVENT_GENERATOR_EVENTGENERATORROOT_H_ +#define MODULES_EVENT_GENERATOR_EVENTGENERATORROOT_H_ + +#include +#include +#include +#include +#include +#include +#include + +#include "../../beans/containers/KinematicRange.h" +#include "EventGeneratorModule.h" + +class TFile; + +class TRandom; + +namespace EPIC { + +/** + * @class EventGeneratorROOT + * + * @brief Event generator based on ROOT library. + * + * This is event generator module based on ROOT library. + */ +class EventGeneratorROOT: public EventGeneratorModule { +public: + static const std::string PARAMETER_NAME_TYPE; ///< Key to set EventGeneratorROOT::m_type. + static const std::string PARAMETER_NAME_ABSTOL; ///< Key to set EventGeneratorROOT::m_absTol. + static const std::string PARAMETER_NAME_RELTOL; ///< Key to set EventGeneratorROOT::m_relTol. + static const std::string PARAMETER_NAME_NBIN; ///< Key to set EventGeneratorROOT::m_nBin. + + static const unsigned int classId; ///< Unique ID to automatically register + /// the class in the registry. + + /** + * Default constructor. + */ + EventGeneratorROOT(const std::string &className); + + /** + * Copy constructor. + */ + EventGeneratorROOT(const EventGeneratorROOT &other); + + /** + * Destructor. + */ + virtual ~EventGeneratorROOT(); + + virtual EventGeneratorROOT *clone() const; + virtual void configure(const ElemUtils::Parameters ¶meters); + + virtual void initialise(const std::vector &kinematicRanges, + const EventGeneratorInterface &service); + virtual std::pair, double> generateEvent(); + virtual std::pair getIntegral(); + +private: + + /** + * Scale value to range. + */ + double scaleValue(const KinematicRange& range, double value) const; + + /** + * Unscale value from range. + */ + double unscaleValue(const KinematicRange& range, double value) const; + + /** + * Get volume based on kinematic ranges. + */ + double getVolume() const; + + ROOT::Math::IntegratorMultiDim *m_pROOT{nullptr}; ///< Pointer to ROOT integrator object. + TRandom *m_pRandom{nullptr}; ///< Pointer to ROOT TRandom object. + TFile* m_rootFile{nullptr}; ///< Pointer to ROOT file. + + const EventGeneratorInterface *m_pEventGeneratorInterface{nullptr}; ///< Pointer used + /// during + /// initialization of ROOT. + std::vector m_kinematicRanges; ///< Kinematic ranges set during initialization. + + std::string m_type{"vegas"}; ///< Type of integration algorithm. + double m_absTol{-1.}; ///< Maximum absolute uncertainty. + double m_relTol{-1.}; ///< Maximum relative uncertainty. + size_t m_nBin{0}; ///< Number of bins in build-up. +}; + +} /* namespace EPIC */ + +#endif /* MODULES_EVENT_GENERATOR_EVENTGENERATORROOT_H_ */ diff --git a/src/modules/event_generator/EventGeneratorROOT.cpp b/src/modules/event_generator/EventGeneratorROOT.cpp new file mode 100644 index 0000000..d0b16e3 --- /dev/null +++ b/src/modules/event_generator/EventGeneratorROOT.cpp @@ -0,0 +1,210 @@ +/* + * EventGeneratorROOT.cpp + * + * Created on: Feb 1, 2024 + * Author: Laurent Forthomme (AGH) + */ + +#include "../../../include/modules/event_generator/EventGeneratorROOT.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../../../include/beans/other/EventGeneratorInterface.h" + +namespace EPIC { + +const std::string EventGeneratorROOT::PARAMETER_NAME_TYPE = "type"; +const std::string EventGeneratorROOT::PARAMETER_NAME_ABSTOL = "absTol"; +const std::string EventGeneratorROOT::PARAMETER_NAME_RELTOL = "relTol"; +const std::string EventGeneratorROOT::PARAMETER_NAME_NBIN = "nBins"; + +const unsigned int EventGeneratorROOT::classId = + PARTONS::BaseObjectRegistry::getInstance()->registerBaseObject( + new EventGeneratorROOT("EventGeneratorROOT")); + +EventGeneratorROOT::EventGeneratorROOT(const std::string &className) : + EventGeneratorModule(className) {} + +EventGeneratorROOT::EventGeneratorROOT(const EventGeneratorROOT &other) : + EventGeneratorModule(other) { + if (other.m_pROOT != nullptr) + warn(__func__, "Not able to copy ROOT object, you need to run " + "initialization for new object"); + + m_type = other.m_type; + m_absTol = other.m_absTol; + m_relTol = other.m_relTol; + m_nBin = other.m_nBin; +} + +EventGeneratorROOT::~EventGeneratorROOT() { + + if (m_pROOT) { + delete m_pROOT; + m_pROOT = nullptr; + } + + if (m_pRandom != nullptr) { + delete m_pRandom; + m_pRandom = nullptr; + } + + if (m_rootFile != nullptr) { + m_rootFile->Close(); + m_rootFile = nullptr; + } +} + +EventGeneratorROOT *EventGeneratorROOT::clone() const { + return new EventGeneratorROOT(*this); +} + +void EventGeneratorROOT::configure(const ElemUtils::Parameters ¶meters) { + + EventGeneratorModule::configure(parameters); + + if (parameters.isAvailable(EventGeneratorROOT::PARAMETER_NAME_TYPE)) { + m_type = parameters.getLastAvailable().toString(); + info(__func__, + ElemUtils::Formatter() << "Parameter " + << EventGeneratorROOT::PARAMETER_NAME_TYPE + << " changed to " << m_type); + } + + if (parameters.isAvailable(EventGeneratorROOT::PARAMETER_NAME_ABSTOL)) { + m_absTol = parameters.getLastAvailable().toDouble(); + info(__func__, + ElemUtils::Formatter() << "Parameter " + << EventGeneratorROOT::PARAMETER_NAME_ABSTOL + << " changed to " << m_absTol); + } + + if (parameters.isAvailable(EventGeneratorROOT::PARAMETER_NAME_RELTOL)) { + m_relTol = parameters.getLastAvailable().toDouble(); + info(__func__, + ElemUtils::Formatter() << "Parameter " + << EventGeneratorROOT::PARAMETER_NAME_RELTOL + << " changed to " << m_relTol); + } + + if (parameters.isAvailable(EventGeneratorROOT::PARAMETER_NAME_NBIN)) { + m_nBin = parameters.getLastAvailable().toUInt(); + info(__func__, + ElemUtils::Formatter() << "Parameter " + << EventGeneratorROOT::PARAMETER_NAME_NBIN + << " changed to " << m_nBin); + } +} + +void EventGeneratorROOT::initialise( + const std::vector &kinematicRanges, + const EventGeneratorInterface &service) { + + //scenario + size_t scenario = 0; + + if (!m_initStatePath.empty()) { + struct stat buffer; //check if file exist + if (!stat(m_initStatePath.c_str(), &buffer) == 0) //no + scenario = 1; + else //yes + scenario = 2; + } + + // random number generator (create) + m_pRandom = new TRandom3(m_seed); + + // kinematic ranges + m_kinematicRanges = kinematicRanges; + + // pointer + m_pEventGeneratorInterface = &service; + + // initialize + if (scenario == 0 || scenario == 1) { + info(__func__, "Creating new ROOT integrator object"); + // parameters + ROOT::Math::IntegratorMultiDim::Type type; + if (m_type == "default") + type = ROOT::Math::IntegratorMultiDim::Type::kDEFAULT; + else if (m_type == "adaptive") + type = ROOT::Math::IntegratorMultiDim::Type::kADAPTIVE; + else if (m_type == "plain") + type = ROOT::Math::IntegratorMultiDim::Type::kPLAIN; + else if (m_type == "miser") + type = ROOT::Math::IntegratorMultiDim::Type::kMISER; + else if (m_type == "vegas") + type = ROOT::Math::IntegratorMultiDim::Type::kVEGAS; + else + throw ElemUtils::CustomException(getClassName(), __func__, + ElemUtils::Formatter() << "Invalid type retrieved: " << m_type); + m_pROOT = new ROOT::Math::IntegratorMultiDim(type, m_absTol, m_relTol, m_nBin); + + if (m_pEventGeneratorInterface == nullptr) + throw ElemUtils::CustomException(getClassName(), __func__, + "Pointer to EventGenerator is null"); + if (m_kinematicRanges.size() != kinematicRanges.size()) + throw ElemUtils::CustomException(getClassName(), __func__, + "Size of vector containing kinematic ranges different than nDim"); + auto density = [&](const double* x) -> double { + std::vector v(kinematicRanges.size()); + for (size_t i = 0; i < kinematicRanges.size(); i++) + v.at(i) = unscaleValue(m_kinematicRanges.at(i), x[i]); + return m_pEventGeneratorInterface->getEventDistribution(v); + }; + m_pROOT->SetFunction(density, kinematicRanges.size() /* number of dimensions */); + } +} + +std::pair, double> EventGeneratorROOT::generateEvent() { + throw ElemUtils::CustomException(getClassName(), __func__, "Module is not yet ready for event generation"); +} + +std::pair EventGeneratorROOT::getIntegral() { + std::vector xmin, xmax; + for (const auto& range : m_kinematicRanges) { + xmin.emplace_back(range.getMin()); + xmax.emplace_back(range.getMax()); + } + if (!m_pROOT) + throw ElemUtils::CustomException(getClassName(), __func__, + "Pointer to ROOT integrator object is null"); + const auto integral = m_pROOT->Integral(xmin.data(), xmax.data()), error = m_pROOT->Error(); + const auto volume = getVolume(); + return std::make_pair(volume * integral, volume * error); +} + +double EventGeneratorROOT::scaleValue(const KinematicRange& range, double value) const { + const auto minScaled = range.getMin(), maxScaled = range.getMax(); + if (value < minScaled || value > maxScaled) + throw ElemUtils::CustomException(getClassName(), __func__, + ElemUtils::Formatter() << "Value outside the limits: " + << value); + return (value - minScaled) / (maxScaled - minScaled); +} + +double EventGeneratorROOT::unscaleValue(const KinematicRange& range, double value) const { + const auto minScaled = range.getMin(), maxScaled = range.getMax(); + if (value < 0. || value > 1.) + throw ElemUtils::CustomException(getClassName(), __func__, + ElemUtils::Formatter() << "Value outside the limits: " << value); + return minScaled + value * (maxScaled - minScaled); +} + +double EventGeneratorROOT::getVolume() const { + double result = 1.; + for (const auto& range : m_kinematicRanges) + result *= range.getMax() - range.getMin(); + return result; +} + +} /* namespace EPIC */ From b7ee6d8c511d4c394a1625937b0b0e72cfcd7c20 Mon Sep 17 00:00:00 2001 From: Laurent Forthomme Date: Fri, 23 Feb 2024 08:38:34 +0100 Subject: [PATCH 2/7] First successful run with a DVCS scenario + ROOT integrator with Miser --- data/examples/epic_scenario_DVCS.xml | 3 + .../event_generator/EventGeneratorROOT.h | 123 +++---- include/services/GeneratorService.h | 2 +- .../event_generator/EventGeneratorROOT.cpp | 317 +++++++----------- 4 files changed, 175 insertions(+), 270 deletions(-) diff --git a/data/examples/epic_scenario_DVCS.xml b/data/examples/epic_scenario_DVCS.xml index 8ca76cf..00643a7 100644 --- a/data/examples/epic_scenario_DVCS.xml +++ b/data/examples/epic_scenario_DVCS.xml @@ -59,6 +59,9 @@ + diff --git a/include/modules/event_generator/EventGeneratorROOT.h b/include/modules/event_generator/EventGeneratorROOT.h index 320c718..525858d 100644 --- a/include/modules/event_generator/EventGeneratorROOT.h +++ b/include/modules/event_generator/EventGeneratorROOT.h @@ -9,9 +9,10 @@ #define MODULES_EVENT_GENERATOR_EVENTGENERATORROOT_H_ #include +#include #include #include -#include + #include #include #include @@ -19,84 +20,44 @@ #include "../../beans/containers/KinematicRange.h" #include "EventGeneratorModule.h" -class TFile; - -class TRandom; - namespace EPIC { - -/** - * @class EventGeneratorROOT - * - * @brief Event generator based on ROOT library. - * - * This is event generator module based on ROOT library. - */ -class EventGeneratorROOT: public EventGeneratorModule { -public: - static const std::string PARAMETER_NAME_TYPE; ///< Key to set EventGeneratorROOT::m_type. - static const std::string PARAMETER_NAME_ABSTOL; ///< Key to set EventGeneratorROOT::m_absTol. - static const std::string PARAMETER_NAME_RELTOL; ///< Key to set EventGeneratorROOT::m_relTol. - static const std::string PARAMETER_NAME_NBIN; ///< Key to set EventGeneratorROOT::m_nBin. - - static const unsigned int classId; ///< Unique ID to automatically register - /// the class in the registry. - - /** - * Default constructor. - */ - EventGeneratorROOT(const std::string &className); - - /** - * Copy constructor. - */ - EventGeneratorROOT(const EventGeneratorROOT &other); - - /** - * Destructor. - */ - virtual ~EventGeneratorROOT(); - - virtual EventGeneratorROOT *clone() const; - virtual void configure(const ElemUtils::Parameters ¶meters); - - virtual void initialise(const std::vector &kinematicRanges, - const EventGeneratorInterface &service); - virtual std::pair, double> generateEvent(); - virtual std::pair getIntegral(); - -private: - - /** - * Scale value to range. - */ - double scaleValue(const KinematicRange& range, double value) const; - - /** - * Unscale value from range. - */ - double unscaleValue(const KinematicRange& range, double value) const; - - /** - * Get volume based on kinematic ranges. - */ - double getVolume() const; - - ROOT::Math::IntegratorMultiDim *m_pROOT{nullptr}; ///< Pointer to ROOT integrator object. - TRandom *m_pRandom{nullptr}; ///< Pointer to ROOT TRandom object. - TFile* m_rootFile{nullptr}; ///< Pointer to ROOT file. - - const EventGeneratorInterface *m_pEventGeneratorInterface{nullptr}; ///< Pointer used - /// during - /// initialization of ROOT. - std::vector m_kinematicRanges; ///< Kinematic ranges set during initialization. - - std::string m_type{"vegas"}; ///< Type of integration algorithm. - double m_absTol{-1.}; ///< Maximum absolute uncertainty. - double m_relTol{-1.}; ///< Maximum relative uncertainty. - size_t m_nBin{0}; ///< Number of bins in build-up. -}; - -} /* namespace EPIC */ - -#endif /* MODULES_EVENT_GENERATOR_EVENTGENERATORROOT_H_ */ + /** + * @class EventGeneratorROOT + * + * @brief Event generator based on ROOT IntegratorMultiDim. + */ + class EventGeneratorROOT : public EventGeneratorModule { + public: + static const std::string PARAMETER_NAME_TYPE; ///< Key to set EventGeneratorROOT::m_type. + static const std::string PARAMETER_NAME_ABSTOL; ///< Key to set EventGeneratorROOT::m_absTol. + static const std::string PARAMETER_NAME_RELTOL; ///< Key to set EventGeneratorROOT::m_relTol. + static const std::string PARAMETER_NAME_NBIN; ///< Key to set EventGeneratorROOT::m_nBin. + + static const unsigned int classId; ///< Unique ID to automatically register the class in the registry. + + EventGeneratorROOT(const std::string &className); ///< Default constructor + EventGeneratorROOT(const EventGeneratorROOT &other); ///< Copy constructor + virtual ~EventGeneratorROOT(); ///< Destructor + + virtual EventGeneratorROOT *clone() const; + virtual void configure(const ElemUtils::Parameters ¶meters); + + virtual void initialise(const std::vector &kinematicRanges, const EventGeneratorInterface &service); + virtual std::pair, double> generateEvent(); + virtual std::pair getIntegral(); + + private: + ROOT::Math::IntegratorMultiDim *m_pROOT{nullptr}; ///< Pointer to ROOT integrator object. + + /// Pointer used during initialization of ROOT. + const EventGeneratorInterface *m_pEventGeneratorInterface{nullptr}; + std::vector m_kinematicRanges; ///< Kinematic ranges set during initialization. + + std::string m_type{"vegas"}; ///< Type of integration algorithm. + double m_absTol{-1.}; ///< Maximum absolute uncertainty. + double m_relTol{-1.}; ///< Maximum relative uncertainty. + size_t m_nBin{0}; ///< Number of bins in build-up. + }; +} // namespace EPIC + +#endif diff --git a/include/services/GeneratorService.h b/include/services/GeneratorService.h index 3203759..1208652 100644 --- a/include/services/GeneratorService.h +++ b/include/services/GeneratorService.h @@ -485,7 +485,7 @@ class GeneratorService: public PARTONS::BaseObject, } if (m_generalConfiguration.getNEvents() == 0) { - throw ElemUtils::CustomException(getClassName(), __func__, + info(__func__, "Number of events to be generated is: 0"); } } diff --git a/src/modules/event_generator/EventGeneratorROOT.cpp b/src/modules/event_generator/EventGeneratorROOT.cpp index d0b16e3..6294f4f 100644 --- a/src/modules/event_generator/EventGeneratorROOT.cpp +++ b/src/modules/event_generator/EventGeneratorROOT.cpp @@ -5,206 +5,147 @@ * Author: Laurent Forthomme (AGH) */ -#include "../../../include/modules/event_generator/EventGeneratorROOT.h" - #include #include #include +#include +#include #include #include -#include #include +#include #include -#include -#include #include "../../../include/beans/other/EventGeneratorInterface.h" +#include "../../../include/modules/event_generator/EventGeneratorROOT.h" namespace EPIC { -const std::string EventGeneratorROOT::PARAMETER_NAME_TYPE = "type"; -const std::string EventGeneratorROOT::PARAMETER_NAME_ABSTOL = "absTol"; -const std::string EventGeneratorROOT::PARAMETER_NAME_RELTOL = "relTol"; -const std::string EventGeneratorROOT::PARAMETER_NAME_NBIN = "nBins"; - -const unsigned int EventGeneratorROOT::classId = - PARTONS::BaseObjectRegistry::getInstance()->registerBaseObject( - new EventGeneratorROOT("EventGeneratorROOT")); - -EventGeneratorROOT::EventGeneratorROOT(const std::string &className) : - EventGeneratorModule(className) {} - -EventGeneratorROOT::EventGeneratorROOT(const EventGeneratorROOT &other) : - EventGeneratorModule(other) { - if (other.m_pROOT != nullptr) - warn(__func__, "Not able to copy ROOT object, you need to run " - "initialization for new object"); - - m_type = other.m_type; - m_absTol = other.m_absTol; - m_relTol = other.m_relTol; - m_nBin = other.m_nBin; -} - -EventGeneratorROOT::~EventGeneratorROOT() { - - if (m_pROOT) { - delete m_pROOT; - m_pROOT = nullptr; - } - - if (m_pRandom != nullptr) { - delete m_pRandom; - m_pRandom = nullptr; - } - - if (m_rootFile != nullptr) { - m_rootFile->Close(); - m_rootFile = nullptr; - } -} - -EventGeneratorROOT *EventGeneratorROOT::clone() const { - return new EventGeneratorROOT(*this); -} - -void EventGeneratorROOT::configure(const ElemUtils::Parameters ¶meters) { - - EventGeneratorModule::configure(parameters); - - if (parameters.isAvailable(EventGeneratorROOT::PARAMETER_NAME_TYPE)) { - m_type = parameters.getLastAvailable().toString(); - info(__func__, - ElemUtils::Formatter() << "Parameter " - << EventGeneratorROOT::PARAMETER_NAME_TYPE - << " changed to " << m_type); - } - - if (parameters.isAvailable(EventGeneratorROOT::PARAMETER_NAME_ABSTOL)) { - m_absTol = parameters.getLastAvailable().toDouble(); - info(__func__, - ElemUtils::Formatter() << "Parameter " - << EventGeneratorROOT::PARAMETER_NAME_ABSTOL - << " changed to " << m_absTol); - } - - if (parameters.isAvailable(EventGeneratorROOT::PARAMETER_NAME_RELTOL)) { - m_relTol = parameters.getLastAvailable().toDouble(); - info(__func__, - ElemUtils::Formatter() << "Parameter " - << EventGeneratorROOT::PARAMETER_NAME_RELTOL - << " changed to " << m_relTol); - } - - if (parameters.isAvailable(EventGeneratorROOT::PARAMETER_NAME_NBIN)) { - m_nBin = parameters.getLastAvailable().toUInt(); - info(__func__, - ElemUtils::Formatter() << "Parameter " - << EventGeneratorROOT::PARAMETER_NAME_NBIN - << " changed to " << m_nBin); - } -} - -void EventGeneratorROOT::initialise( - const std::vector &kinematicRanges, - const EventGeneratorInterface &service) { - - //scenario - size_t scenario = 0; - - if (!m_initStatePath.empty()) { - struct stat buffer; //check if file exist - if (!stat(m_initStatePath.c_str(), &buffer) == 0) //no - scenario = 1; - else //yes - scenario = 2; - } - - // random number generator (create) - m_pRandom = new TRandom3(m_seed); - - // kinematic ranges - m_kinematicRanges = kinematicRanges; - - // pointer - m_pEventGeneratorInterface = &service; - - // initialize - if (scenario == 0 || scenario == 1) { - info(__func__, "Creating new ROOT integrator object"); - // parameters - ROOT::Math::IntegratorMultiDim::Type type; - if (m_type == "default") - type = ROOT::Math::IntegratorMultiDim::Type::kDEFAULT; - else if (m_type == "adaptive") - type = ROOT::Math::IntegratorMultiDim::Type::kADAPTIVE; - else if (m_type == "plain") - type = ROOT::Math::IntegratorMultiDim::Type::kPLAIN; - else if (m_type == "miser") - type = ROOT::Math::IntegratorMultiDim::Type::kMISER; - else if (m_type == "vegas") - type = ROOT::Math::IntegratorMultiDim::Type::kVEGAS; - else - throw ElemUtils::CustomException(getClassName(), __func__, - ElemUtils::Formatter() << "Invalid type retrieved: " << m_type); - m_pROOT = new ROOT::Math::IntegratorMultiDim(type, m_absTol, m_relTol, m_nBin); - - if (m_pEventGeneratorInterface == nullptr) - throw ElemUtils::CustomException(getClassName(), __func__, - "Pointer to EventGenerator is null"); - if (m_kinematicRanges.size() != kinematicRanges.size()) - throw ElemUtils::CustomException(getClassName(), __func__, - "Size of vector containing kinematic ranges different than nDim"); - auto density = [&](const double* x) -> double { - std::vector v(kinematicRanges.size()); - for (size_t i = 0; i < kinematicRanges.size(); i++) - v.at(i) = unscaleValue(m_kinematicRanges.at(i), x[i]); - return m_pEventGeneratorInterface->getEventDistribution(v); - }; - m_pROOT->SetFunction(density, kinematicRanges.size() /* number of dimensions */); - } -} - -std::pair, double> EventGeneratorROOT::generateEvent() { - throw ElemUtils::CustomException(getClassName(), __func__, "Module is not yet ready for event generation"); -} - -std::pair EventGeneratorROOT::getIntegral() { - std::vector xmin, xmax; - for (const auto& range : m_kinematicRanges) { - xmin.emplace_back(range.getMin()); - xmax.emplace_back(range.getMax()); + const std::string EventGeneratorROOT::PARAMETER_NAME_TYPE = "type"; + const std::string EventGeneratorROOT::PARAMETER_NAME_ABSTOL = "absTol"; + const std::string EventGeneratorROOT::PARAMETER_NAME_RELTOL = "relTol"; + const std::string EventGeneratorROOT::PARAMETER_NAME_NBIN = "nBins"; + + const unsigned int EventGeneratorROOT::classId = + PARTONS::BaseObjectRegistry::getInstance()->registerBaseObject(new EventGeneratorROOT("EventGeneratorROOT")); + + EventGeneratorROOT::EventGeneratorROOT(const std::string &className) : EventGeneratorModule(className) {} + + EventGeneratorROOT::EventGeneratorROOT(const EventGeneratorROOT &other) : EventGeneratorModule(other) { + if (other.m_pROOT != nullptr) + warn(__func__, + "Not able to copy ROOT object, you need to run " + "initialization for new object"); + + m_type = other.m_type; + m_absTol = other.m_absTol; + m_relTol = other.m_relTol; + m_nBin = other.m_nBin; + } + + EventGeneratorROOT::~EventGeneratorROOT() { + if (m_pROOT) + delete m_pROOT; + } + + EventGeneratorROOT *EventGeneratorROOT::clone() const { return new EventGeneratorROOT(*this); } + + void EventGeneratorROOT::configure(const ElemUtils::Parameters ¶meters) { + EventGeneratorModule::configure(parameters); + + if (parameters.isAvailable(EventGeneratorROOT::PARAMETER_NAME_TYPE)) { + m_type = parameters.getLastAvailable().getString(); + info(__func__, + ElemUtils::Formatter() << "Parameter " << EventGeneratorROOT::PARAMETER_NAME_TYPE << " changed to " + << m_type); + } + + if (parameters.isAvailable(EventGeneratorROOT::PARAMETER_NAME_ABSTOL)) { + m_absTol = parameters.getLastAvailable().toDouble(); + info(__func__, + ElemUtils::Formatter() << "Parameter " << EventGeneratorROOT::PARAMETER_NAME_ABSTOL << " changed to " + << m_absTol); + } + + if (parameters.isAvailable(EventGeneratorROOT::PARAMETER_NAME_RELTOL)) { + m_relTol = parameters.getLastAvailable().toDouble(); + info(__func__, + ElemUtils::Formatter() << "Parameter " << EventGeneratorROOT::PARAMETER_NAME_RELTOL << " changed to " + << m_relTol); + } + + if (parameters.isAvailable(EventGeneratorROOT::PARAMETER_NAME_NBIN)) { + m_nBin = parameters.getLastAvailable().toUInt(); + info(__func__, + ElemUtils::Formatter() << "Parameter " << EventGeneratorROOT::PARAMETER_NAME_NBIN << " changed to " + << m_nBin); + } + } + + void EventGeneratorROOT::initialise(const std::vector &kinematicRanges, + const EventGeneratorInterface &service) { + //scenario + size_t scenario = 0; + + if (!m_initStatePath.empty()) { + struct stat buffer; + if (!stat(m_initStatePath.c_str(), &buffer) == 0) // file does not exist + scenario = 1; + else // file exists + scenario = 2; + } + + m_kinematicRanges = kinematicRanges; + m_pEventGeneratorInterface = &service; + + // initialize + if (scenario == 0 || scenario == 1) { + info(__func__, "Creating new ROOT integrator object"); + // parameters + ROOT::Math::IntegratorMultiDim::Type type; + if (m_type == "default") + type = ROOT::Math::IntegratorMultiDim::Type::kDEFAULT; + else if (m_type == "adaptive") + type = ROOT::Math::IntegratorMultiDim::Type::kADAPTIVE; + else if (m_type == "plain") + type = ROOT::Math::IntegratorMultiDim::Type::kPLAIN; + else if (m_type == "miser") + type = ROOT::Math::IntegratorMultiDim::Type::kMISER; + else if (m_type == "vegas") + type = ROOT::Math::IntegratorMultiDim::Type::kVEGAS; + else + throw ElemUtils::CustomException( + getClassName(), __func__, ElemUtils::Formatter() << "Invalid type retrieved: " << m_type); + m_pROOT = new ROOT::Math::IntegratorMultiDim(type, m_absTol, m_relTol, m_nBin); + + if (m_pEventGeneratorInterface == nullptr) + throw ElemUtils::CustomException(getClassName(), __func__, "Pointer to EventGenerator is null"); + if (m_kinematicRanges.size() != kinematicRanges.size()) + throw ElemUtils::CustomException( + getClassName(), __func__, "Size of vector containing kinematic ranges different than nDim"); + static const auto density = [&](const double *x) -> double { + std::vector v(kinematicRanges.size()); + for (size_t i = 0; i < kinematicRanges.size(); i++) + v.at(i) = x[i]; + return m_pEventGeneratorInterface->getEventDistribution(v); + }; + m_pROOT->SetFunction(density, kinematicRanges.size() /* number of dimensions */); + } + } + + std::pair, double> EventGeneratorROOT::generateEvent() { + throw ElemUtils::CustomException(getClassName(), __func__, "Module is not yet ready for event generation"); + } + + std::pair EventGeneratorROOT::getIntegral() { + std::vector xmin, xmax; + for (const auto &range : m_kinematicRanges) { + xmin.emplace_back(range.getMin()); + xmax.emplace_back(range.getMax()); + } + if (!m_pROOT) + throw ElemUtils::CustomException(getClassName(), __func__, "Pointer to ROOT integrator object is null"); + return std::make_pair(m_pROOT->Integral(xmin.data(), xmax.data()), m_pROOT->Error()); } - if (!m_pROOT) - throw ElemUtils::CustomException(getClassName(), __func__, - "Pointer to ROOT integrator object is null"); - const auto integral = m_pROOT->Integral(xmin.data(), xmax.data()), error = m_pROOT->Error(); - const auto volume = getVolume(); - return std::make_pair(volume * integral, volume * error); -} - -double EventGeneratorROOT::scaleValue(const KinematicRange& range, double value) const { - const auto minScaled = range.getMin(), maxScaled = range.getMax(); - if (value < minScaled || value > maxScaled) - throw ElemUtils::CustomException(getClassName(), __func__, - ElemUtils::Formatter() << "Value outside the limits: " - << value); - return (value - minScaled) / (maxScaled - minScaled); -} - -double EventGeneratorROOT::unscaleValue(const KinematicRange& range, double value) const { - const auto minScaled = range.getMin(), maxScaled = range.getMax(); - if (value < 0. || value > 1.) - throw ElemUtils::CustomException(getClassName(), __func__, - ElemUtils::Formatter() << "Value outside the limits: " << value); - return minScaled + value * (maxScaled - minScaled); -} - -double EventGeneratorROOT::getVolume() const { - double result = 1.; - for (const auto& range : m_kinematicRanges) - result *= range.getMax() - range.getMin(); - return result; -} } /* namespace EPIC */ From 1f3401145c1bfc07b81bd8de722f98e23f07984e Mon Sep 17 00:00:00 2001 From: Laurent Forthomme Date: Fri, 23 Feb 2024 13:43:36 +0100 Subject: [PATCH 3/7] Fix on uncertainty retrieval --- src/modules/event_generator/EventGeneratorROOT.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/modules/event_generator/EventGeneratorROOT.cpp b/src/modules/event_generator/EventGeneratorROOT.cpp index 6294f4f..07f97a3 100644 --- a/src/modules/event_generator/EventGeneratorROOT.cpp +++ b/src/modules/event_generator/EventGeneratorROOT.cpp @@ -145,7 +145,8 @@ namespace EPIC { } if (!m_pROOT) throw ElemUtils::CustomException(getClassName(), __func__, "Pointer to ROOT integrator object is null"); - return std::make_pair(m_pROOT->Integral(xmin.data(), xmax.data()), m_pROOT->Error()); + const auto integral = m_pROOT->Integral(xmin.data(), xmax.data()), uncertainty = m_pROOT->Error(); + return std::make_pair(integral, uncertainty); } } /* namespace EPIC */ From 0065dec90234686e816d29f98a297d693c4c4b2c Mon Sep 17 00:00:00 2001 From: Laurent Forthomme Date: Fri, 23 Feb 2024 14:25:08 +0100 Subject: [PATCH 4/7] Setting and and printout of integrator parameters --- src/modules/event_generator/EventGeneratorROOT.cpp | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/modules/event_generator/EventGeneratorROOT.cpp b/src/modules/event_generator/EventGeneratorROOT.cpp index 07f97a3..0f568e6 100644 --- a/src/modules/event_generator/EventGeneratorROOT.cpp +++ b/src/modules/event_generator/EventGeneratorROOT.cpp @@ -117,6 +117,14 @@ namespace EPIC { throw ElemUtils::CustomException( getClassName(), __func__, ElemUtils::Formatter() << "Invalid type retrieved: " << m_type); m_pROOT = new ROOT::Math::IntegratorMultiDim(type, m_absTol, m_relTol, m_nBin); + m_pROOT->Options().SetAbsTolerance(m_absTol); + m_pROOT->Options().SetRelTolerance(m_relTol); + m_pROOT->Options().SetWKSize(m_nBin); + { + std::ostringstream os; + m_pROOT->Options().Print(os); + info(__func__, ElemUtils::Formatter() << "List of options for ROOT integrator\n" << os.str()); + } if (m_pEventGeneratorInterface == nullptr) throw ElemUtils::CustomException(getClassName(), __func__, "Pointer to EventGenerator is null"); From a35715e05fcb22a0e763f44c894c595859bc11a7 Mon Sep 17 00:00:00 2001 From: Laurent Forthomme Date: Wed, 13 Mar 2024 17:28:33 +0100 Subject: [PATCH 5/7] Perform integration before (any) event generation --- .../event_generator/EventGeneratorROOT.h | 2 ++ .../event_generator/EventGeneratorROOT.cpp | 21 +++++++++---------- 2 files changed, 12 insertions(+), 11 deletions(-) diff --git a/include/modules/event_generator/EventGeneratorROOT.h b/include/modules/event_generator/EventGeneratorROOT.h index 525858d..e9ede2a 100644 --- a/include/modules/event_generator/EventGeneratorROOT.h +++ b/include/modules/event_generator/EventGeneratorROOT.h @@ -57,6 +57,8 @@ namespace EPIC { double m_absTol{-1.}; ///< Maximum absolute uncertainty. double m_relTol{-1.}; ///< Maximum relative uncertainty. size_t m_nBin{0}; ///< Number of bins in build-up. + + std::pair m_integral{0., 0.}; ///< Integrated cross section + uncertainty }; } // namespace EPIC diff --git a/src/modules/event_generator/EventGeneratorROOT.cpp b/src/modules/event_generator/EventGeneratorROOT.cpp index 0f568e6..d17e270 100644 --- a/src/modules/event_generator/EventGeneratorROOT.cpp +++ b/src/modules/event_generator/EventGeneratorROOT.cpp @@ -138,6 +138,15 @@ namespace EPIC { return m_pEventGeneratorInterface->getEventDistribution(v); }; m_pROOT->SetFunction(density, kinematicRanges.size() /* number of dimensions */); + std::vector xmin, xmax; + for (const auto &range : m_kinematicRanges) { + xmin.emplace_back(range.getMin()); + xmax.emplace_back(range.getMax()); + } + if (!m_pROOT) + throw ElemUtils::CustomException(getClassName(), __func__, "Pointer to ROOT integrator object is null"); + const auto integral = m_pROOT->Integral(xmin.data(), xmax.data()), uncertainty = m_pROOT->Error(); + m_integral = std::make_pair(integral, uncertainty); } } @@ -145,16 +154,6 @@ namespace EPIC { throw ElemUtils::CustomException(getClassName(), __func__, "Module is not yet ready for event generation"); } - std::pair EventGeneratorROOT::getIntegral() { - std::vector xmin, xmax; - for (const auto &range : m_kinematicRanges) { - xmin.emplace_back(range.getMin()); - xmax.emplace_back(range.getMax()); - } - if (!m_pROOT) - throw ElemUtils::CustomException(getClassName(), __func__, "Pointer to ROOT integrator object is null"); - const auto integral = m_pROOT->Integral(xmin.data(), xmax.data()), uncertainty = m_pROOT->Error(); - return std::make_pair(integral, uncertainty); - } + std::pair EventGeneratorROOT::getIntegral() { return m_integral; } } /* namespace EPIC */ From 6c54b98404f2cdf01fdcb9b9189e8a57568e0db2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Pawe=C5=82=20Sznajder?= <88622027+pawelsznajder@users.noreply.github.com> Date: Wed, 13 Mar 2024 17:49:18 +0100 Subject: [PATCH 6/7] Update epic_scenario_DVCS.xml clean --- data/examples/epic_scenario_DVCS.xml | 3 --- 1 file changed, 3 deletions(-) diff --git a/data/examples/epic_scenario_DVCS.xml b/data/examples/epic_scenario_DVCS.xml index 00643a7..8ca76cf 100644 --- a/data/examples/epic_scenario_DVCS.xml +++ b/data/examples/epic_scenario_DVCS.xml @@ -59,9 +59,6 @@ - From 433fa0c3c6a3bc226ed37007dc9cb52dc514d5f3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Pawe=C5=82=20Sznajder?= <88622027+pawelsznajder@users.noreply.github.com> Date: Wed, 13 Mar 2024 17:50:52 +0100 Subject: [PATCH 7/7] Update GeneratorService.h change info to warn --- include/services/GeneratorService.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/services/GeneratorService.h b/include/services/GeneratorService.h index 1208652..0cbdd35 100644 --- a/include/services/GeneratorService.h +++ b/include/services/GeneratorService.h @@ -485,7 +485,7 @@ class GeneratorService: public PARTONS::BaseObject, } if (m_generalConfiguration.getNEvents() == 0) { - info(__func__, + warn(__func__, "Number of events to be generated is: 0"); } }