From 7c8c7beae453284278dc03d9627f620956943325 Mon Sep 17 00:00:00 2001 From: AndreaRigoni Date: Thu, 19 Mar 2026 14:57:26 +0000 Subject: [PATCH] add ecomug --- src/HEP/Geant/EcoMug.hh | 1203 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 1203 insertions(+) create mode 100644 src/HEP/Geant/EcoMug.hh diff --git a/src/HEP/Geant/EcoMug.hh b/src/HEP/Geant/EcoMug.hh new file mode 100644 index 0000000..1271b1f --- /dev/null +++ b/src/HEP/Geant/EcoMug.hh @@ -0,0 +1,1203 @@ +///////////////////////////////////////////////////////////////////////////////////// +// EcoMug: Efficient COsmic MUon Generator // +// Copyright (C) 2023 Davide Pagano // +// // +// EcoMug is based on the following work: // +// D. Pagano, G. Bonomi, A. Donzella, A. Zenoni, G. Zumerle, N. Zurlo, // +// "EcoMug: an Efficient COsmic MUon Generator for cosmic-ray muons applications", // +// doi:10.1016/j.nima.2021.165732 // +// // +// 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 // +// the Free Software Foundation, either version 3 of the License, or // +// (at your option) any later version. // +// // +// This program is distributed in the hope that it will be useful, // +// but WITHOUT ANY WARRANTY; without even the implied warranty of // +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // +// GNU General Public License for more details. // +// // +// You should have received a copy of the GNU General Public License // +// along with this program. If not, see . // +///////////////////////////////////////////////////////////////////////////////////// + +#ifndef EcoMug_H +#define EcoMug_H + +#include +#include +#include +#include +#include +#include +#include + +#define ECOMUG_VERSION "2.1" + +#ifndef M_PI +# define M_PI_NOT_DEFINED +# define M_PI 3.14159265358979323846 +#endif + +namespace EMUnits { + // Default units: + // meter (m) + // second (s) + // Giga electron Volt (GeV) + // radian (rad) + + // Lengths and areas + static const double m = 1.; + static const double cm = 1.e-2*m; + static const double mm = 1.e-3*m; + static const double km = 1000.*m; + static const double mm2 = mm*mm; + static const double cm2 = cm*cm; + static const double m2 = m*m; + static const double km2 = km*km; + + // Angles + static const double rad = 1.; + static const double mrad = 1.e-3*rad; + static const double deg = (M_PI/180.0)*rad; + + // Time + static const double s = 1.; + static const double ms = 1.e-3*s; + static const double us = 1.e-6*s; + static const double ns = 1.e-9*s; + static const double min = 60.*s; + static const double hour = 60.*min; + static const double day = 24.*hour; + static const double hertz = 1./s; + + // Energy/momentum + static const double GeV = 1.; + static const double MeV = 1.e-3*GeV; + static const double keV = 1.e-3*MeV; + static const double TeV = 1.e+6*MeV; + static const double eV = 1.e-6*MeV; +}; +/////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////// + + +//! Class for the logging system +class EMLog { +public: + enum TLogLevel {ERROR, WARNING, INFO, DEBUG}; + + enum TMsgType {EcoMug, EMRandom, EMMultiGen, EMMaximization, UNKNOWN}; + + EMLog() {}; + virtual ~EMLog() { + os << std::endl; + fprintf(stderr, "%s", os.str().c_str()); + fflush(stderr); + }; + + std::ostringstream& Get(TLogLevel level = INFO); + void Get(TLogLevel level, const std::string& msg, TMsgType type = EcoMug) { + std::cout << "[EcoMug v" << ECOMUG_VERSION << "] ["; + std::cout << ToString(level); + if (type != UNKNOWN) std::cout << " in " << ToString(type); + std::cout << "]: "; + std::cout << std::string(level > DEBUG ? level - DEBUG : 0, '\t'); + std::cout << msg << "\n"; + }; + +public: + static std::string ToString(TLogLevel level) { + static const char* const buffer[] = {"ERROR", "WARNING", "INFO", "DEBUG"}; + return buffer[level]; + }; + + static std::string ToString(TMsgType type) { + static const char* const buffer[] = {"EcoMug", "EMRandom", "EMMultiGen", "EMMaximization", "UNKNOWN"}; + return buffer[type]; + }; + + static TLogLevel FromString(const std::string& level) { + if (level == "DEBUG") + return DEBUG; + if (level == "INFO") + return INFO; + if (level == "WARNING") + return WARNING; + if (level == "ERROR") + return ERROR; + EMLog().Get(WARNING) << "Unknown logging level '" << level << "'. Using INFO level as default."; + return INFO; + }; + + static EMLog::TLogLevel ReportingLevel; +private: + EMLog(const EMLog&); + EMLog& operator =(const EMLog&); + std::ostringstream os; +}; + +#define EMLogger(level, msg, type) \ +if (level > EMLog::ReportingLevel) ; \ +else EMLog().Get(level, msg, type) + +//EMLog::TLogLevel EMLog::ReportingLevel = WARNING; +/////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////// + + +//! Fast generation of random numbers +//! This class is based on the xoroshiro128+ generator. +//! https://prng.di.unimi.it/ +class EMRandom { +public: + EMRandom() { + std::random_device rd; + std::mt19937 gen(rd()); + std::uniform_int_distribution dis(0, std::numeric_limits::max()); + s[0] = dis(gen); + s[1] = dis(gen); + }; + + void SetSeed(std::uint64_t seed) { + s[0] = seed; + s[1] = seed; + }; + + double GenerateRandomDouble() { + std::uint64_t x = next(); + return to_double(x); + }; + + double GenerateRandomDouble(double x1, double x2) { + return (x2-x1)*GenerateRandomDouble()+x1; + }; + + int64_t rotl(const std::uint64_t x, int k) { + return (x << k) | (x >> (64 - k)); + }; + + std::uint64_t next() { + const std::uint64_t s0 = s[0]; + std::uint64_t s1 = s[1]; + const std::uint64_t result = s0 + s1; + s1 ^= s0; + s[0] = rotl(s0, 55) ^ s1 ^ (s1 << 14); + s[1] = rotl(s1, 36); + return result; + }; + + double to_double(std::uint64_t x) { + union U { + std::uint64_t i; + double d; + }; + U u = { UINT64_C(0x3FF) << 52 | x >> 12 }; + return u.d - 1.0; + }; + + std::uint64_t s[2]; +}; +/////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////// + + +//! Class for maximization based on "Whale Optimization Algorithm" +//! doi:10.1016/j.advengsoft.2016.01.008 +class EMMaximization { +private: + EMRandom mRandom; + std::size_t mPopSize; + std::size_t mNIter; + int mGenMethod; // 0 = sky, 1 = cylinder, 2 = hspere + double m_a; + double m_a2; + std::vector > mRanges; + std::vector > mPopulation; + double mBestCost; + std::vector mBestSolution; + std::function mFunc; + +public: + EMMaximization(const EMRandom& random, int genMethod) : mRandom(random), mPopSize(200), + mNIter(500), mGenMethod(genMethod), m_a(0.), m_a2(0.), mBestCost(-1.) { + mFunc = &DefaultJ; + }; + /////////////////////////////////////////////////////////////// + + static double DefaultJ(double p, double theta) { + double n = std::max(0.1, 2.856-0.655*log(p)); + return 1600*pow(p, 0.279)*pow(cos(theta), n); + } + /////////////////////////////////////////////////////////////// + + void SetParameters(double minP, double maxP, double minTheta, double maxTheta) { + mRanges.push_back({minP, maxP}); + mRanges.push_back({minTheta, maxTheta}); + } + /////////////////////////////////////////////////////////////// + + void SetParameters(double minP, double maxP, double minTheta, double maxTheta, double minPhi, double maxPhi) { + mRanges.push_back({minP, maxP}); + mRanges.push_back({minTheta, maxTheta}); + mRanges.push_back({minPhi, maxPhi}); + mRanges.push_back({0, M_PI/2.}); + } + /////////////////////////////////////////////////////////////// + + void SetFunction(std::function func) { + mFunc = func; + } + /////////////////////////////////////////////////////////////// + + double SkyFunc(double p, double theta) { + return mFunc(p, theta)*cos(theta)*sin(theta); + } + /////////////////////////////////////////////////////////////// + + double CylFunc(double p, double theta) { + return mFunc(p, theta)*pow(sin(theta), 2); + } + /////////////////////////////////////////////////////////////// + + double HSFunc(double p, double theta, double phi, double theta0) { + return mFunc(p, theta)*(sin(theta0)*sin(theta)*cos(phi) + cos(theta0)*cos(theta))*sin(theta); + } + /////////////////////////////////////////////////////////////// + + double Evaluate(std::vector &v) { + if (mGenMethod == 0) { + return SkyFunc(v[0], v[1]); + } else if (mGenMethod == 1) { + return CylFunc(v[0], v[1]); + } else { + return HSFunc(v[0], v[1], v[2], v[3]); + } + return -1; + } + /////////////////////////////////////////////////////////////// + + void Evaluate() { + double value; + for (std::size_t i = 0; i < mPopSize; ++i) { + value = Evaluate(mPopulation[i]); + if (value > mBestCost) { + mBestCost = value; + mBestSolution = mPopulation[i]; + } + } + } + /////////////////////////////////////////////////////////////// + + void Init() { + std::size_t dim = mRanges.size(); + mPopulation.resize(mPopSize); + for (std::size_t i = 0; i < mPopSize; ++i) { + mPopulation[i].resize(dim); + for (std::size_t j = 0; j < dim; ++j) { + mPopulation[i][j] = mRandom.GenerateRandomDouble(mRanges[j][0], mRanges[j][1]); + } + } + } + /////////////////////////////////////////////////////////////// + + void UpdateParameters(std::size_t t) { + m_a = 2. - t*(2./mNIter); + m_a2 = -1. + t*((-1.)/mNIter); + } + /////////////////////////////////////////////////////////////// + + void Move() { + double r1, r2, A, C, b, l, rw, p, D_tmp, D_best, distance; + std::vector tmp; + for (std::size_t i = 0; i < mPopulation.size(); ++i) { + r1 = mRandom.GenerateRandomDouble(); + r2 = mRandom.GenerateRandomDouble(); + A = 2*m_a*r1-m_a; + C = 2*r2; + b = 1.; + l = (m_a2-1)*mRandom.GenerateRandomDouble()+1; + p = mRandom.GenerateRandomDouble(); + + for (std::size_t j = 0; j < mPopulation[0].size(); ++j) { + if (p < 0.5) { + if (fabs(A) >= 1) { + rw = floor(mRandom.GenerateRandomDouble()*mPopulation.size()); + tmp = mPopulation[rw]; + D_tmp = fabs(C*tmp[j] - mPopulation[i][j]); + mPopulation[i][j] = tmp[j] - A*D_tmp; + } else { + D_best = fabs(C*mBestSolution[j] - mPopulation[i][j]); + mPopulation[i][j] = mBestSolution[j]-A*D_best; + } + } else { + distance = fabs(mBestSolution[j] - mPopulation[i][j]); + mPopulation[i][j] = distance*exp(b*l)*cos(l*2*M_PI) + mBestSolution[j]; + } + if (mPopulation[i][j] < mRanges[j][0]) mPopulation[i][j] = mRanges[j][0]; + if (mPopulation[i][j] > mRanges[j][1]) mPopulation[i][j] = mRanges[j][1]; + } + } + } + /////////////////////////////////////////////////////////////// + + double Maximize() { + Init(); + Evaluate(); + for (std::size_t iter = 1; iter < mNIter; ++iter) { + UpdateParameters(iter); + Move(); + Evaluate(); + } + return mBestCost; + } +}; +/////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////// + + +//! Class for the generation of cosmic muons +class EcoMug { + friend class EMRandom; + +public: + /// Possible generation methods + enum EMGeometry { + Sky, ///< generation from a plane (flat sky) + Cylinder, ///< generation from a cylinder + HSphere ///< generation from a half-sphere + }; + +private: + EMGeometry mGenMethod; + std::array mGenerationPosition; + double mGenerationTheta; + double mGenerationPhi; + double mGenerationMomentum; + double mMinimumMomentum; + double mMaximumMomentum; + double mMinimumTheta; + double mMaximumTheta; + double mMinimumPhi; + double mMaximumPhi; + int mCharge; + double mHorizontalRate; + double mCylinderMinPositionPhi; + double mCylinderMaxPositionPhi; + double mHSphereMinPositionPhi; + double mHSphereMaxPositionPhi; + double mHSphereMinPositionTheta; + double mHSphereMaxPositionTheta; + double mHSphereCosMinPositionTheta; + double mHSphereCosMaxPositionTheta; + double mJPrime; + double mN; + double mRandAccRej; + double mPhi0; + double mTheta0; + bool mAccepted; + std::array mSkySize; + std::array mSkyCenterPosition; + double mCylinderHeight; + double mCylinderRadius; + std::array mCylinderCenterPosition; + double mHSphereRadius; + double mMaxFuncSkyCylinder; + std::array mHSphereCenterPosition; + bool mCustomJ; + EMRandom mRandom; + std::default_random_engine mEngineC; + std::discrete_distribution mDiscDistC; + std::array mMaxJ; + std::array mMaxCustomJ; + std::function mJ; + +public: + // Default constructor + EcoMug() : mGenMethod(Sky), + mGenerationPosition({{0., 0., 0.}}), mGenerationTheta(0.), mGenerationPhi(0.), + mGenerationMomentum(0.), mMinimumMomentum(0.01), mMaximumMomentum(1000.), + mMinimumTheta(0.), mMaximumTheta(M_PI/2.), mMinimumPhi(0.), mMaximumPhi(2.*M_PI), + mCharge(1), mHorizontalRate(129*EMUnits::hertz/EMUnits::m2), mCylinderMinPositionPhi(0.), mCylinderMaxPositionPhi(2.*M_PI), + mHSphereMinPositionPhi(0.), mHSphereMaxPositionPhi(2.*M_PI), mHSphereMinPositionTheta(0.), + mHSphereMaxPositionTheta(M_PI/2.), mHSphereCosMinPositionTheta(1.), mHSphereCosMaxPositionTheta(0.), + mJPrime(0.), mN(0.), mRandAccRej(0.), mPhi0(0.), mTheta0(0.), mAccepted(false), + mSkySize({{0., 0.}}), mSkyCenterPosition({{0., 0., 0.}}), mCylinderHeight(0.), + mCylinderRadius(0.), mCylinderCenterPosition({{0., 0., 0.}}), mHSphereRadius(0.), + mMaxFuncSkyCylinder(5.3176), mHSphereCenterPosition({{0., 0., 0.}}), mCustomJ(false), + mEngineC(std::random_device{}()) { + mDiscDistC = std::discrete_distribution({128, 100}); + mMaxJ = {-1., -1., -1.}; + mMaxCustomJ = {-1., -1., -1.}; + }; + + // Copy constructor + EcoMug(const EcoMug& t) { + mGenMethod = t.mGenMethod; + mGenerationPosition = t.mGenerationPosition; + mGenerationTheta = t.mGenerationTheta; + mGenerationPhi = t.mGenerationPhi; + mGenerationMomentum = t.mGenerationMomentum; + mMinimumMomentum = t.mMinimumMomentum; + mMaximumMomentum = t.mMaximumMomentum; + mMinimumTheta = t.mMinimumTheta; + mMaximumTheta = t.mMaximumTheta; + mMinimumPhi = t.mMinimumPhi; + mMaximumPhi = t.mMaximumPhi; + mCharge = t.mCharge; + mHorizontalRate = t.mHorizontalRate; + mCylinderMinPositionPhi = t.mCylinderMinPositionPhi; + mCylinderMaxPositionPhi = t.mCylinderMaxPositionPhi; + mHSphereMinPositionPhi = t.mHSphereMinPositionPhi; + mHSphereMaxPositionPhi = t.mHSphereMaxPositionPhi; + mHSphereMinPositionTheta = t.mHSphereMinPositionTheta; + mHSphereMaxPositionTheta = t.mHSphereMaxPositionTheta; + mHSphereCosMinPositionTheta = t.mHSphereCosMinPositionTheta; + mHSphereCosMaxPositionTheta = t.mHSphereCosMaxPositionTheta; + mJPrime = t.mJPrime; + mN = t.mN; + mRandAccRej = t.mRandAccRej; + mPhi0 = t.mPhi0; + mTheta0 = t.mTheta0; + mAccepted = t.mAccepted; + mSkySize = t.mSkySize; + mSkyCenterPosition = t.mSkyCenterPosition; + mCylinderHeight = t.mCylinderHeight; + mCylinderRadius = t.mCylinderRadius; + mCylinderCenterPosition = t.mCylinderCenterPosition; + mHSphereRadius = t.mHSphereRadius; + mMaxFuncSkyCylinder = t.mMaxFuncSkyCylinder; + mHSphereCenterPosition = t.mHSphereCenterPosition; + mCustomJ = t.mCustomJ; + mRandom = t.mRandom; + mEngineC = t.mEngineC; + mDiscDistC = t.mDiscDistC; + mMaxJ = t.mMaxJ; + mMaxCustomJ = t.mMaxCustomJ; + mJ = t.mJ; + }; + + + /////////////////////////////////////////////////////////////// + // Methods to access the parameters of the generated muon + /////////////////////////////////////////////////////////////// + /// Get the generation position + const std::array& GetGenerationPosition() const { + return mGenerationPosition; + }; + /// Get the generation momentum + double GetGenerationMomentum() const { + return mGenerationMomentum; + }; + /// Get the generation momentum + void GetGenerationMomentum(std::array& momentum) const { + momentum = { + mGenerationMomentum*sin(mGenerationTheta)*cos(mGenerationPhi), + mGenerationMomentum*sin(mGenerationTheta)*sin(mGenerationPhi), + mGenerationMomentum*cos(mGenerationTheta) + }; + }; + /// Get the generation theta + double GetGenerationTheta() const { + return mGenerationTheta; + }; + /// Get the generation phi + double GetGenerationPhi() const { + return mGenerationPhi; + }; + /// Get charge + int GetCharge() const { + return mCharge; + }; + /////////////////////////////////////////////////////////////// + + + /////////////////////////////////////////////////////////////// + // Methods for the geometry of the generation + /////////////////////////////////////////////////////////////// + /// Set generation from sky + void SetUseSky() { + mGenMethod = Sky; + }; + /// Set cylindrical generation + void SetUseCylinder() { + mGenMethod = Cylinder; + }; + /// Set half-sphere generation + void SetUseHSphere() { + mGenMethod = HSphere; + }; + /// Set the generation method (Sky, Cylinder or HSphere) + void SetGenerationMethod(EMGeometry genM) { + mGenMethod = genM; + }; + + /// Get the generation method (Sky, Cylinder or HSphere) + EMGeometry GetGenerationMethod() const { + return mGenMethod; + }; + /////////////////////////////////////////////////////////////// + + + /////////////////////////////////////////////////////////////// + // Common methods to all geometries + /////////////////////////////////////////////////////////////// + /// Set the differential flux J. Accepted functions are like + /// double J(double momentum, double theta) + /// momentum has to be in GeV/c and theta in radians + void SetDifferentialFlux(std::function J) { + mJ = J; + mCustomJ = true; + }; + /// Set the seed for the internal PRNG (if 0 a random seed is used) + void SetSeed(std::uint64_t seed) { + if (seed > 0) mRandom.SetSeed(seed); + }; + /// Set minimum generation Momentum + void SetMinimumMomentum(double momentum) { + mMinimumMomentum = momentum; + }; + /// Set maximum generation Momentum + void SetMaximumMomentum(double momentum) { + mMaximumMomentum = momentum; + }; + /// Set minimum generation Theta + void SetMinimumTheta(double theta) { + mMinimumTheta = theta; + }; + /// Set maximum generation Theta + void SetMaximumTheta(double theta) { + mMaximumTheta = theta; + }; + /// Set minimum generation Phi + void SetMinimumPhi(double phi) { + mMinimumPhi = phi; + }; + /// Set maximum generation Phi + void SetMaximumPhi(double phi) { + mMaximumPhi = phi; + }; + /// Set the rate of cosmic ray muons per square unit are through + /// a horizontal surface. Default value is 129 Hz/m^2. + void SetHorizontalRate(double rate) { + mHorizontalRate = rate; + }; + + /// Get minimum generation Momentum + double GetMinimumMomentum() const { + return mMinimumMomentum; + }; + /// Get maximum generation Momentum + double GetMaximumMomentum() const { + return mMaximumMomentum; + }; + /// Get minimum generation Theta + double GetMinimumTheta() const { + return mMinimumTheta; + }; + /// Get maximum generation Theta + double GetMaximumTheta() const { + return mMaximumTheta; + }; + /// Get minimum generation Phi + double GetMinimumPhi() const { + return mMinimumPhi; + }; + /// Get maximum generation Phi + double GetMaximumPhi() const { + return mMaximumPhi; + }; + /// Get horizontal rate of cosmic muons + double GetHorizontalRate() const { + return mHorizontalRate; + }; + /// Get the generation surface area + double GetGenSurfaceArea() const { + double area = 0.; + if (mGenMethod == Sky) { + area = mSkySize[0]*mSkySize[1]; + } else if (mGenMethod == Cylinder) { + // A = \Delta\phi r h + area = (mCylinderMaxPositionPhi-mCylinderMinPositionPhi)*mCylinderRadius*mCylinderHeight; + } else { + // A = \Delta\phi r^2\left(\cos\theta_{min} - \cos\theta_{max}\right) + area = (mHSphereMaxPositionPhi-mHSphereMinPositionPhi)*pow(mHSphereRadius, 2)*(cos(mHSphereMinPositionTheta) - cos(mHSphereMaxPositionTheta)); + } + return area; + }; + /// Get the average rate of cosmic muons as well the error on it + /// The optional parameter npoints defines the number + /// of points to be used in the MC integration of the J' + void GetAverageGenRateAndError(double &rate, double &error, int npoints = 1e7) { + // 129.0827 is the integral of the J'prime for the sky in + // the full range of theta, phi and up to 3 TeV in energy. + // For custom J it is the user who should account for the correction. + double k = mHorizontalRate/129.0827; + if (mGenMethod == Sky) { + if (mCustomJ) { + MCJprimeCustomSkyIntegration(rate, error, npoints); + return; + } else MCJprimeSkyIntegration(rate, error, npoints); + } else if (mGenMethod == Cylinder) { + if (mCustomJ) { + MCJprimeCustomCylinderIntegration(rate, error, npoints); + return; + } else MCJprimeCylinderIntegration(rate, error, npoints); + } else { + if (mCustomJ) { + MCJprimeCustomHSphereIntegration(rate, error, npoints); + return; + } else MCJprimeHSphereIntegration(rate, error, npoints); + } + rate *= k; + error *= k; + }; + /// Get the average rate of cosmic muons + /// The optional parameter npoints defines the number + /// of points to be used in the MC integration of the J' + double GetAverageGenRate(int npoints = 1e7) { + double rate, error; + GetAverageGenRateAndError(rate, error, npoints); + return rate; + }; + /// Get the estimated corresponding to the provided statistics + double GetEstimatedTime(int nmuons) { + if (mCustomJ) return 0.; + return (nmuons/(GetGenSurfaceArea()/EMUnits::m2))/(GetAverageGenRate()/EMUnits::hertz*EMUnits::m2); + }; + /////////////////////////////////////////////////////////////// + + + /////////////////////////////////////////////////////////////// + // Methods for the plane-based generation + /////////////////////////////////////////////////////////////// + /// Set sky size + void SetSkySize(const std::array& size) { + mSkySize = size; + }; + /// Set sky center position + void SetSkyCenterPosition(const std::array& position) { + mSkyCenterPosition = position; + }; + double GetSkySize(unsigned int index) { + return mSkySize[index]; + }; + /////////////////////////////////////////////////////////////// + + + /////////////////////////////////////////////////////////////// + // Methods for the cylinder-based generation + /////////////////////////////////////////////////////////////// + /// Set cylinder radius + void SetCylinderRadius(double radius) { + mCylinderRadius = radius; + }; + /// Set cylinder height + void SetCylinderHeight(double height) { + mCylinderHeight = height; + }; + /// Set cylinder center position + void SetCylinderCenterPosition(const std::array& position) { + mCylinderCenterPosition = position; + }; + void SetCylinderMinPositionPhi(double phi) { + mCylinderMinPositionPhi = phi; + }; + void SetCylinderMaxPositionPhi(double phi) { + mCylinderMaxPositionPhi = phi; + }; + /// Get cylinder radius + double GetCylinderRadius() const { + return mCylinderRadius; + }; + /// Get cylinder height + double GetCylinderHeight() const { + return mCylinderHeight; + }; + /// Get cylinder center position + const std::array& GetCylinderCenterPosition() const { + return mCylinderCenterPosition; + }; + /////////////////////////////////////////////////////////////// + + + /////////////////////////////////////////////////////////////// + // Methods for the half sphere-based generation + /////////////////////////////////////////////////////////////// + /// Set half-sphere radius + void SetHSphereRadius(double radius) { + mHSphereRadius = radius; + }; + /// Set half-sphere center position + void SetHSphereCenterPosition(const std::array& position) { + mHSphereCenterPosition = position; + }; + void SetHSphereMinPositionPhi(double phi) { + mHSphereMinPositionPhi = phi; + }; + void SetHSphereMaxPositionPhi(double phi) { + mHSphereMaxPositionPhi = phi; + }; + void SetHSphereMinPositionTheta(double theta) { + mHSphereMinPositionTheta = theta; + mHSphereCosMinPositionTheta = cos(mHSphereMinPositionTheta); + }; + void SetHSphereMaxPositionTheta(double theta) { + mHSphereMaxPositionTheta = theta; + mHSphereCosMaxPositionTheta = cos(mHSphereMaxPositionTheta); + }; + /// Get half-sphere radius + double GetHSphereRadius() const { + return mHSphereRadius; + }; + /// Get half-sphere center position + const std::array& GetHSphereCenterPosition() const { + return mHSphereCenterPosition; + }; + /////////////////////////////////////////////////////////////// + + +private: + double F1Cumulative(double x) { + return 1. - 8.534790171171021/pow(x + 2.68, 87./40.); + }; + + double F1Inverse(double x) { + return (2.68 - 2.68*pow(1. - x, 40./87.))/pow(1. - x, 40./87.); + }; + + double maxSkyJFunc() { + return 1600*pow(mMaximumMomentum, 0.279)*pow(cos(0.76158), 1.1)*sin(0.76158); + }; + + double maxCylJFunc() { + return 1600*pow(mMaximumMomentum, 0.279)*pow(cos(1.35081), 0.1)*pow(sin(1.35081), 2); + }; + + double maxHSJFunc() { + return 1600*pow(mMaximumMomentum, 0.279)*pow(cos(1.26452), 0.1)*(sin(1.26452)*sin(1.26452)+cos(1.26452)*cos(1.26452))*sin(1.26452); + }; + + double GenerateMomentumF1() { + double z = mRandom.GenerateRandomDouble(F1Cumulative(mMinimumMomentum), F1Cumulative(mMaximumMomentum)); + return F1Inverse(z); + }; + + void GeneratePositionSky() { + mGenerationPosition[0] = mRandom.GenerateRandomDouble(mSkyCenterPosition[0]-mSkySize[0]/2., mSkyCenterPosition[0]+mSkySize[0]/2.); + mGenerationPosition[1] = mRandom.GenerateRandomDouble(mSkyCenterPosition[1]-mSkySize[1]/2., mSkyCenterPosition[1]+mSkySize[1]/2.); + mGenerationPosition[2] = mSkyCenterPosition[2]; + }; + + void GeneratePositionCylinder() { + mPhi0 = mRandom.GenerateRandomDouble(mCylinderMinPositionPhi, mCylinderMaxPositionPhi); + mGenerationPosition[0] = mCylinderCenterPosition[0] + mCylinderRadius*cos(mPhi0); + mGenerationPosition[1] = mCylinderCenterPosition[1] + mCylinderRadius*sin(mPhi0); + mGenerationPosition[2] = mRandom.GenerateRandomDouble(mCylinderCenterPosition[2]-mCylinderHeight/2., mCylinderCenterPosition[2]+mCylinderHeight/2.); + }; + + void ComputeMaximumCustomJ() { + EMMaximization maximizer(mRandom, mGenMethod); + maximizer.SetFunction(mJ); + if (mGenMethod == 0 || mGenMethod == 1) { + maximizer.SetParameters(mMinimumMomentum, mMaximumMomentum, mMinimumTheta, mMaximumTheta); + } else { + maximizer.SetParameters(mMinimumMomentum, mMaximumMomentum, mMinimumTheta, mMaximumTheta, mMinimumPhi, mMaximumPhi); + } + mMaxCustomJ[mGenMethod] = maximizer.Maximize(); + }; + + void ComputeMaximum() { + EMMaximization maximizer(mRandom, mGenMethod); + if (mGenMethod == 0 || mGenMethod == 1) { + maximizer.SetParameters(mMinimumMomentum, mMaximumMomentum, mMinimumTheta, mMaximumTheta); + } else { + maximizer.SetParameters(mMinimumMomentum, mMaximumMomentum, mMinimumTheta, mMaximumTheta, mMinimumPhi, mMaximumPhi); + } + mMaxJ[mGenMethod] = maximizer.Maximize(); + }; + + void MCJprimeCustomSkyIntegration(double &rate, double &error, int npoints) { + double I = 0., I2 = 0., value = 0.; + for (auto i = 0; i < npoints; ++i) { + mGenerationTheta = mRandom.GenerateRandomDouble(mMinimumTheta, mMaximumTheta); + mGenerationMomentum = mRandom.GenerateRandomDouble(mMinimumMomentum, mMaximumMomentum); + value = mJ(mGenerationMomentum, mGenerationTheta)*cos(mGenerationTheta)*sin(mGenerationTheta); + I += value; + I2 += pow(value, 2); + } + double V = (mMaximumMomentum-mMinimumMomentum)*(mMaximumTheta-mMinimumTheta)*(mMaximumPhi-mMinimumPhi); + double expected = I/npoints; + double expectedSquare = I2/npoints; + rate = V*I/npoints; + error = V*pow((expectedSquare-pow(expected,2))/(npoints-1), 0.5); + }; + + void MCJprimeCustomCylinderIntegration(double &rate, double &error, int npoints) { + double I = 0., I2 = 0., value = 0.; + for (auto i = 0; i < npoints; ++i) { + mGenerationTheta = mRandom.GenerateRandomDouble(mMinimumTheta, mMaximumTheta); + mGenerationPhi = mRandom.GenerateRandomDouble(mMinimumPhi, mMaximumPhi); + mGenerationMomentum = mRandom.GenerateRandomDouble(mMinimumMomentum, mMaximumMomentum); + value = mJ(mGenerationMomentum, mGenerationTheta)*pow(sin(mGenerationTheta), 2)*cos(mGenerationPhi); + if (value < 0) value = 0; + I += value; + I2 += pow(value, 2); + } + double V = (mMaximumMomentum-mMinimumMomentum)*(mMaximumTheta-mMinimumTheta)*(mMaximumPhi-mMinimumPhi); + double expected = I/npoints; + double expectedSquare = I2/npoints; + rate = V*I/npoints; + error = V*pow((expectedSquare-pow(expected,2))/(npoints-1), 0.5); + }; + + void MCJprimeCustomHSphereIntegration(double &rate, double &error, int npoints) { + double I = 0., I2 = 0., value = 0.; + for (auto i = 0; i < npoints; ++i) { + mTheta0 = mRandom.GenerateRandomDouble(mHSphereMinPositionTheta, mHSphereMaxPositionTheta); + mGenerationTheta = mRandom.GenerateRandomDouble(mMinimumTheta, mMaximumTheta); + mGenerationPhi = mRandom.GenerateRandomDouble(mMinimumPhi, mMaximumPhi); + mGenerationMomentum = mRandom.GenerateRandomDouble(mMinimumMomentum, mMaximumMomentum); + value = mJ(mGenerationMomentum, mGenerationTheta)*(sin(mTheta0)*sin(mGenerationTheta)*sin(mGenerationTheta)*cos(mGenerationPhi)+cos(mTheta0)*cos(mGenerationTheta)*sin(mGenerationTheta))*sin(mTheta0); + if (value < 0) value = 0; + I += value; + I2 += pow(value, 2); + } + double V = (mMaximumMomentum-mMinimumMomentum)*(mMaximumTheta-mMinimumTheta)*(mMaximumPhi-mMinimumPhi)*(mHSphereMaxPositionTheta-mHSphereMinPositionTheta); + double expected = I/npoints; + double expectedSquare = I2/npoints; + rate = V*I/npoints; + error = V*pow((expectedSquare-pow(expected,2))/(npoints-1), 0.5); + }; + + void MCJprimeSkyIntegration(double &rate, double &error, int npoints) { + double I = 0., I2 = 0., value = 0.; + for (auto i = 0; i < npoints; ++i) { + mGenerationTheta = mRandom.GenerateRandomDouble(mMinimumTheta, mMaximumTheta); + mGenerationMomentum = mRandom.GenerateRandomDouble(mMinimumMomentum, mMaximumMomentum); + mN = 2.856-0.655*log(mGenerationMomentum); + if (mN < 0.1) mN = 0.1; + value = 1600*pow(mGenerationMomentum+2.68, -3.175)*pow(mGenerationMomentum, 0.279)*pow(cos(mGenerationTheta), mN+1)*sin(mGenerationTheta); + I += value; + I2 += pow(value, 2); + } + double V = (mMaximumMomentum-mMinimumMomentum)*(mMaximumTheta-mMinimumTheta)*(mMaximumPhi-mMinimumPhi); + double expected = I/npoints; + double expectedSquare = I2/npoints; + rate = V*I/npoints; + error = V*pow((expectedSquare-pow(expected,2))/(npoints-1), 0.5); + }; + + void MCJprimeCylinderIntegration(double &rate, double &error, int npoints) { + double I = 0., I2 = 0., value = 0.; + for (auto i = 0; i < npoints; ++i) { + mGenerationTheta = mRandom.GenerateRandomDouble(mMinimumTheta, mMaximumTheta); + mGenerationPhi = mRandom.GenerateRandomDouble(mMinimumPhi, mMaximumPhi); + mGenerationMomentum = mRandom.GenerateRandomDouble(mMinimumMomentum, mMaximumMomentum); + mN = 2.856-0.655*log(mGenerationMomentum); + if (mN < 0.1) mN = 0.1; + value = 1600*pow(mGenerationMomentum+2.68, -3.175)*pow(mGenerationMomentum, 0.279)*pow(cos(mGenerationTheta), mN)*pow(sin(mGenerationTheta), 2)*cos(mGenerationPhi); + if (value < 0) value = 0; + I += value; + I2 += pow(value, 2); + } + double V = (mMaximumMomentum-mMinimumMomentum)*(mMaximumTheta-mMinimumTheta)*(mMaximumPhi-mMinimumPhi); + double expected = I/npoints; + double expectedSquare = I2/npoints; + rate = V*I/npoints; + error = V*pow((expectedSquare-pow(expected,2))/(npoints-1), 0.5); + }; + + void MCJprimeHSphereIntegration(double &rate, double &error, int npoints) { + double I = 0., I2 = 0., value = 0.; + for (auto i = 0; i < npoints; ++i) { + mTheta0 = mRandom.GenerateRandomDouble(mHSphereMinPositionTheta, mHSphereMaxPositionTheta); + mGenerationTheta = mRandom.GenerateRandomDouble(mMinimumTheta, mMaximumTheta); + mGenerationPhi = mRandom.GenerateRandomDouble(mMinimumPhi, mMaximumPhi); + mGenerationMomentum = mRandom.GenerateRandomDouble(mMinimumMomentum, mMaximumMomentum); + mN = 2.856-0.655*log(mGenerationMomentum); + if (mN < 0.1) mN = 0.1; + value = 1600*pow(mGenerationMomentum+2.68, -3.175)*pow(mGenerationMomentum, 0.279)*pow(cos(mGenerationTheta), mN)*(sin(mTheta0)*sin(mGenerationTheta)*sin(mGenerationTheta)*cos(mGenerationPhi)+cos(mTheta0)*cos(mGenerationTheta)*sin(mGenerationTheta))*sin(mTheta0); + if (value < 0) value = 0; + I += value; + I2 += pow(value, 2); + } + double V = (mMaximumMomentum-mMinimumMomentum)*(mMaximumTheta-mMinimumTheta)*(mMaximumPhi-mMinimumPhi)*(mHSphereMaxPositionTheta-mHSphereMinPositionTheta); + double expected = I/npoints; + double expectedSquare = I2/npoints; + rate = V*I/npoints; + error = V*pow((expectedSquare-pow(expected,2))/(npoints-1), 0.5); + }; + +public: + /////////////////////////////////////////////////////////////// + /// Generate a cosmic muon from the pre-defined J + /////////////////////////////////////////////////////////////// + void Generate() { + mAccepted = false; + + if (mMaxJ[mGenMethod] < 0) ComputeMaximum(); + + // Sky or cylinder generation + if (mGenMethod == Sky || mGenMethod == Cylinder) { + // Generation of the momentum and theta angle + while (!mAccepted) { + mRandAccRej = mRandom.GenerateRandomDouble(); + mGenerationTheta = mRandom.GenerateRandomDouble(mMinimumTheta, mMaximumTheta); + mGenerationMomentum = GenerateMomentumF1(); + mN = 2.856-0.655*log(mGenerationMomentum); + if (mN < 0.1) mN = 0.1; + + if (mGenMethod == Sky) { + mJPrime = 1600*pow(mGenerationMomentum, 0.279)*pow(cos(mGenerationTheta), mN+1)*sin(mGenerationTheta); + if (mMaxJ[mGenMethod]*mRandAccRej < mJPrime) mAccepted = true; + } + + if(mGenMethod == Cylinder) { + mJPrime = 1600*pow(mGenerationMomentum, 0.279)*pow(cos(mGenerationTheta), mN)*pow(sin(mGenerationTheta), 2); + if (mMaxJ[mGenMethod]*mRandAccRej < mJPrime) mAccepted = true; + } + } + mGenerationTheta = M_PI - mGenerationTheta; + + // Generation of the position and phi angle + if (mGenMethod == Sky) { + GeneratePositionSky(); + mGenerationPhi = mRandom.GenerateRandomDouble(mMinimumPhi, mMaximumPhi); + } + if (mGenMethod == Cylinder) { + mAccepted = false; + GeneratePositionCylinder(); + while (!mAccepted) { + mRandAccRej = mRandom.GenerateRandomDouble(); + mGenerationPhi = mRandom.GenerateRandomDouble(mMinimumPhi, mMaximumPhi); + if (mRandAccRej < fabs(cos(mGenerationPhi))) mAccepted = true; + } + mGenerationPhi = mGenerationPhi + mPhi0; + if (mGenerationPhi >= 2.*M_PI) mGenerationPhi -= 2.*M_PI; + + // Check if the muon is inward + if (sin(mGenerationTheta)*cos(mGenerationPhi)*mGenerationPosition[0] + sin(mGenerationTheta)*sin(mGenerationPhi)*mGenerationPosition[1] > 0) Generate(); + } + } + + // Half-sphere generation + if (mGenMethod == HSphere) { + // Generation point on the half-sphere + mPhi0 = mRandom.GenerateRandomDouble(mHSphereMinPositionPhi, mHSphereMaxPositionPhi); + while (!mAccepted) { + mRandAccRej = mRandom.GenerateRandomDouble(); + mTheta0 = acos(mRandom.GenerateRandomDouble(mHSphereCosMaxPositionTheta, mHSphereCosMinPositionTheta)); + mGenerationTheta = mRandom.GenerateRandomDouble(mMinimumTheta, mMaximumTheta); + mGenerationPhi = mRandom.GenerateRandomDouble(mMinimumPhi, mMaximumPhi); + mGenerationMomentum = GenerateMomentumF1(); + mN = 2.856-0.655*log(mGenerationMomentum); + if (mN < 0.1) mN = 0.1; + + mJPrime = 1600*pow(mGenerationMomentum, 0.279)*pow(cos(mGenerationTheta), mN)*(sin(mGenerationTheta)*sin(mTheta0)*cos(mGenerationPhi)+cos(mGenerationTheta)*cos(mTheta0))*sin(mGenerationTheta); + if (mJPrime > 0 && mMaxJ[mGenMethod]*mRandAccRej < mJPrime) mAccepted = true; + } + + mGenerationPosition[0] = mHSphereRadius*sin(mTheta0)*cos(mPhi0) + mHSphereCenterPosition[0]; + mGenerationPosition[1] = mHSphereRadius*sin(mTheta0)*sin(mPhi0) + mHSphereCenterPosition[1]; + mGenerationPosition[2] = mHSphereRadius*cos(mTheta0) + mHSphereCenterPosition[2]; + + mGenerationTheta = M_PI - mGenerationTheta; + mGenerationPhi = mGenerationPhi + mPhi0; + if (mGenerationPhi >= 2*M_PI) mGenerationPhi -= 2*M_PI; + + mGenerationPhi += M_PI; + if (mGenerationPhi >= 2*M_PI) mGenerationPhi -= 2*M_PI; + } + + // Generate the charge + mCharge = (mDiscDistC(mEngineC) == 0) ? 1 : -1; + }; + /////////////////////////////////////////////////////////////// + + + /////////////////////////////////////////////////////////////// + /// Generate a cosmic muon for the user-defined J + /////////////////////////////////////////////////////////////// + void GenerateFromCustomJ() { + mAccepted = false; + + std::cout << mMaxCustomJ[mGenMethod] << std::endl; + + if (mMaxCustomJ[mGenMethod] < 0) ComputeMaximumCustomJ(); + + // Sky or cylinder generation + if (mGenMethod == Sky || mGenMethod == Cylinder) { + // Generation of the momentum and theta angle + while (!mAccepted) { + mRandAccRej = mRandom.GenerateRandomDouble(); + mGenerationTheta = mRandom.GenerateRandomDouble(mMinimumTheta, mMaximumTheta); + mGenerationMomentum = mRandom.GenerateRandomDouble(mMinimumMomentum, mMaximumMomentum); + + if (mGenMethod == Sky) { + mJPrime = mJ(mGenerationMomentum, mGenerationTheta)*cos(mGenerationTheta)*sin(mGenerationTheta); + if (mMaxCustomJ[mGenMethod]*mRandAccRej < mJPrime) mAccepted = true; + } + + if(mGenMethod == Cylinder) { + mJPrime = mJ(mGenerationMomentum, mGenerationTheta)*pow(sin(mGenerationTheta), 2)*cos(mGenerationPhi); + if (mMaxCustomJ[mGenMethod]*mRandAccRej < mJPrime) mAccepted = true; + } + } + mGenerationTheta = M_PI - mGenerationTheta; + + // Generation of the position and phi angle + if (mGenMethod == Sky) { + GeneratePositionSky(); + mGenerationPhi = mRandom.GenerateRandomDouble(mMinimumPhi, mMaximumPhi); + } + if (mGenMethod == Cylinder) { + mAccepted = false; + GeneratePositionCylinder(); + while (!mAccepted) { + mRandAccRej = mRandom.GenerateRandomDouble(); + mGenerationPhi = mRandom.GenerateRandomDouble(mMinimumPhi, mMaximumPhi); + if (mRandAccRej < fabs(cos(mGenerationPhi))) mAccepted = true; + } + mGenerationPhi = mGenerationPhi + mPhi0; + if (mGenerationPhi >= 2.*M_PI) mGenerationPhi -= 2.*M_PI; + + // Check if the muon is inward + if (sin(mGenerationTheta)*cos(mGenerationPhi)*mGenerationPosition[0] + sin(mGenerationTheta)*sin(mGenerationPhi)*mGenerationPosition[1] > 0) Generate(); + } + } + + // Half-sphere generation + if (mGenMethod == HSphere) { + // Generation point on the half-sphere + mPhi0 = mRandom.GenerateRandomDouble(mHSphereMinPositionPhi, mHSphereMaxPositionPhi); + while (!mAccepted) { + mRandAccRej = mRandom.GenerateRandomDouble(); + mTheta0 = acos(mRandom.GenerateRandomDouble(mHSphereCosMaxPositionTheta, mHSphereCosMinPositionTheta)); + mGenerationTheta = mRandom.GenerateRandomDouble(mMinimumTheta, mMaximumTheta); + mGenerationPhi = mRandom.GenerateRandomDouble(mMinimumPhi, mMaximumPhi); + mGenerationMomentum = mRandom.GenerateRandomDouble(mMinimumMomentum, mMaximumMomentum); + + mJPrime = mJ(mGenerationMomentum, mGenerationTheta)*(sin(mTheta0)*sin(mGenerationTheta)*cos(mGenerationPhi) + cos(mTheta0)*cos(mGenerationTheta))*sin(mGenerationTheta); + if (mMaxCustomJ[mGenMethod]*mRandAccRej < mJPrime) mAccepted = true; + } + + mGenerationPosition[0] = mHSphereRadius*sin(mTheta0)*cos(mPhi0) + mHSphereCenterPosition[0]; + mGenerationPosition[1] = mHSphereRadius*sin(mTheta0)*sin(mPhi0) + mHSphereCenterPosition[1]; + mGenerationPosition[2] = mHSphereRadius*cos(mTheta0) + mHSphereCenterPosition[2]; + + mGenerationTheta = M_PI - mGenerationTheta; + mGenerationPhi = mGenerationPhi + mPhi0; + if (mGenerationPhi >= 2*M_PI) mGenerationPhi -= 2*M_PI; + + mGenerationPhi += M_PI; + if (mGenerationPhi >= 2*M_PI) mGenerationPhi -= 2*M_PI; + } + + // Generate the charge + mCharge = (mDiscDistC(mEngineC) == 0) ? 1 : -1; + }; + /////////////////////////////////////////////////////////////// + +}; +/////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////// + + +//! Class to handle the generation from multiple distributions, +//! for example to take into account backgound sources. +class EMMultiGen { +public: +EMMultiGen(const EcoMug& signal, const std::vector& backgrounds) : + mIndex(-1), mSigInstance(signal), mBckInstances{backgrounds}, mLimits(backgrounds.size()+2), mWeights(backgrounds.size()+1, 1.), + mPID(backgrounds.size()+1), mRd(std::random_device{}()) { + for (std::size_t i = 0; i < mLimits.size(); ++i) mLimits[i] = i; + mDd = std::piecewise_constant_distribution<>(mLimits.begin(), mLimits.end(), mWeights.begin()); +}; + +/// Set the weights for all EcoMug background instance. The number of elements +/// must be equal to the background instances defined +void SetBckWeights(const std::vector& weights) { + if (mBckInstances.size() != weights.size()) { + EMLogger(EMLog::ERROR, "Expected " + std::to_string(mBckInstances.size()) + " weights, but " + std::to_string(weights.size()) + " were provided. Setting them to 1.", EMLog::EMMultiGen); + std::fill(mWeights.begin(), mWeights.end(), 1); + } else { + for (std::size_t i = 0; i < weights.size(); ++i) mWeights[i+1] = weights[i]; + } + mDd = std::piecewise_constant_distribution<>(mLimits.begin(), mLimits.end(), mWeights.begin()); +}; + +/// Set the PID for all background instances. +void SetBckPID(const std::vector& values) { + if (mBckInstances.size() != values.size()) { + EMLogger(EMLog::ERROR, "Expected " + std::to_string(mBckInstances.size()) + " PID, but " + std::to_string(values.size()) + " were provided. Setting them to 0.", EMLog::EMMultiGen); + std::fill(mPID.begin(), mPID.end(), 0); + return; + } + for (std::size_t i = 0; i < values.size(); ++i) mPID[i+1] = values[i]; +}; + +/// Get the generation position +const std::array& GetGenerationPosition() const { + return mBckInstances[mIndex].GetGenerationPosition(); +}; + +/// Get the generation momentum +double GetGenerationMomentum() const { + return mBckInstances[mIndex].GetGenerationMomentum(); +}; + +/// Get the generation momentum +void GetGenerationMomentum(std::array& momentum) const { + mBckInstances[mIndex].GetGenerationMomentum(momentum); +}; + +/// Get the generation theta +double GetGenerationTheta() const { + return mBckInstances[mIndex].GetGenerationTheta(); +}; + +/// Get the generation phi +double GetGenerationPhi() const { + return mBckInstances[mIndex].GetGenerationPhi(); +}; + +/// Get PID +int GetPID() const { + // muon case + if (mPID[mIndex] == 0) { + if (mSigInstance.GetCharge() < 0) return 13; + else return -13; + } + return mPID[mIndex]; +}; + +void Generate() { + mIndex = (int) mDd(mRd); + if (mIndex == 0) mSigInstance.Generate(); + else mBckInstances[mIndex-1].Generate(); +}; + +private: +int mIndex; +EcoMug mSigInstance; +std::vector mBckInstances; +std::vector mLimits; +std::vector mWeights; +std::vector mPID; +std::default_random_engine mRd; +std::piecewise_constant_distribution<> mDd; +}; +/////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////// + +#ifdef ECOMUG_VERSION +#undef ECOMUG_VERSION +#endif + +#ifdef M_PI_NOT_DEFINED +#undef M_PI +#undef M_PI_NOT_DEFINED +#endif + +#endif \ No newline at end of file