diff --git a/DDCore/src/ShapeTags.h b/DDCore/include/DD4hep/ShapeTags.h similarity index 100% rename from DDCore/src/ShapeTags.h rename to DDCore/include/DD4hep/ShapeTags.h diff --git a/DDCore/src/ShapeUtilities.cpp b/DDCore/src/ShapeUtilities.cpp index b3fdd84d57895ad9d2c9d731e621807497a0af7d..61a5aed5ffd25ce3f55d31e0bda8a409f193986a 100644 --- a/DDCore/src/ShapeUtilities.cpp +++ b/DDCore/src/ShapeUtilities.cpp @@ -14,12 +14,12 @@ #define _USE_MATH_DEFINES // Framework include files -#include "DD4hep/Shapes.h" -#include "DD4hep/MatrixHelpers.h" -#include "DD4hep/DD4hepUnits.h" -#include "DD4hep/Printout.h" -#include "DD4hep/detail/ShapesInterna.h" -#include "ShapeTags.h" +#include <DD4hep/Shapes.h> +#include <DD4hep/MatrixHelpers.h> +#include <DD4hep/DD4hepUnits.h> +#include <DD4hep/ShapeTags.h> +#include <DD4hep/Printout.h> +#include <DD4hep/detail/ShapesInterna.h> // C/C++ include files #include <stdexcept> @@ -27,11 +27,11 @@ #include <sstream> // ROOT includes -#include "TClass.h" -#include "TGeoMatrix.h" -#include "TGeoBoolNode.h" -#include "TGeoScaledShape.h" -#include "TGeoCompositeShape.h" +#include <TClass.h> +#include <TGeoMatrix.h> +#include <TGeoBoolNode.h> +#include <TGeoScaledShape.h> +#include <TGeoCompositeShape.h> using namespace std; diff --git a/DDCore/src/Shapes.cpp b/DDCore/src/Shapes.cpp index 8b47177c2cade8365ab878cfb0347e6e73b9ce3f..3dd2fc58317eda7c921dcb08ce4624b2a4637461 100644 --- a/DDCore/src/Shapes.cpp +++ b/DDCore/src/Shapes.cpp @@ -14,12 +14,12 @@ #define _USE_MATH_DEFINES // Framework include files -#include "DD4hep/Detector.h" -#include "DD4hep/MatrixHelpers.h" -#include "DD4hep/DD4hepUnits.h" -#include "DD4hep/Printout.h" -#include "DD4hep/detail/ShapesInterna.h" -#include "ShapeTags.h" +#include <DD4hep/Detector.h> +#include <DD4hep/MatrixHelpers.h> +#include <DD4hep/DD4hepUnits.h> +#include <DD4hep/ShapeTags.h> +#include <DD4hep/Printout.h> +#include <DD4hep/detail/ShapesInterna.h> // C/C++ include files #include <stdexcept> @@ -27,11 +27,11 @@ #include <sstream> // ROOT includes -#include "TClass.h" -#include "TGeoMatrix.h" -#include "TGeoBoolNode.h" -#include "TGeoScaledShape.h" -#include "TGeoCompositeShape.h" +#include <TClass.h> +#include <TGeoMatrix.h> +#include <TGeoBoolNode.h> +#include <TGeoScaledShape.h> +#include <TGeoCompositeShape.h> using namespace std; using namespace dd4hep; diff --git a/DDCore/src/gdml/DetElementCreator.cpp b/DDCore/src/gdml/DetElementCreator.cpp index ffdbe144fd0e3a2c6ddead84d9bcd23f4fba2cd4..598217e4a770db0cb2dbac827b20a06289b47d43 100644 --- a/DDCore/src/gdml/DetElementCreator.cpp +++ b/DDCore/src/gdml/DetElementCreator.cpp @@ -96,7 +96,7 @@ namespace dd4hep { int sd_lvl, PrintLevel p); /// Default destructor - virtual ~DetElementCreator(); + virtual ~DetElementCreator() throw(); /// Callback to output PlacedVolume information of an single Placement virtual int operator()(PlacedVolume pv, int level); /// Callback to output PlacedVolume information of an entire Placement @@ -140,7 +140,8 @@ DetElementCreator::DetElementCreator(Detector& desc, } /// Default destructor -DetElementCreator::~DetElementCreator() { +DetElementCreator::~DetElementCreator() throw() +{ Count total; stringstream str, id_str; const char* pref = detector_volume_match.c_str(); diff --git a/DDCore/src/plugins/Compact2Objects.cpp b/DDCore/src/plugins/Compact2Objects.cpp index c38958bb1abd811795354fcc83c7e431da425826..9522469fdb05923c5352aa72302cc37c72c3fbf5 100644 --- a/DDCore/src/plugins/Compact2Objects.cpp +++ b/DDCore/src/plugins/Compact2Objects.cpp @@ -18,34 +18,34 @@ //========================================================================== // // Framework includes -#include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/DetectorTools.h" -#include "DD4hep/MatrixHelpers.h" -#include "DD4hep/PropertyTable.h" -#include "DD4hep/OpticalSurfaces.h" -#include "DD4hep/OpticalSurfaceManager.h" -#include "DD4hep/IDDescriptor.h" -#include "DD4hep/DD4hepUnits.h" -#include "DD4hep/FieldTypes.h" -#include "DD4hep/Printout.h" -#include "DD4hep/Plugins.h" -#include "DD4hep/detail/SegmentationsInterna.h" -#include "DD4hep/detail/DetectorInterna.h" -#include "DD4hep/detail/ObjectsInterna.h" - -#include "XML/DocumentHandler.h" -#include "XML/Utilities.h" +#include <DD4hep/DetFactoryHelper.h> +#include <DD4hep/DetectorTools.h> +#include <DD4hep/MatrixHelpers.h> +#include <DD4hep/PropertyTable.h> +#include <DD4hep/OpticalSurfaces.h> +#include <DD4hep/OpticalSurfaceManager.h> +#include <DD4hep/IDDescriptor.h> +#include <DD4hep/DD4hepUnits.h> +#include <DD4hep/FieldTypes.h> +#include <DD4hep/Printout.h> +#include <DD4hep/Plugins.h> +#include <DD4hep/detail/SegmentationsInterna.h> +#include <DD4hep/detail/DetectorInterna.h> +#include <DD4hep/detail/ObjectsInterna.h> + +#include <XML/DocumentHandler.h> +#include <XML/Utilities.h> // Root/TGeo include files -#include "TGeoManager.h" -#include "TGeoMaterial.h" +#include <TGeoManager.h> +#include <TGeoMaterial.h> #if ROOT_VERSION_CODE >= ROOT_VERSION(6,12,0) -#include "TGeoPhysicalConstants.h" +#include <TGeoPhysicalConstants.h> #endif #if ROOT_VERSION_CODE >= ROOT_VERSION(6,17,0) -#include "TGDMLMatrix.h" +#include <TGDMLMatrix.h> #endif -#include "TMath.h" +#include <TMath.h> // C/C++ include files #include <climits> diff --git a/DDCore/src/plugins/ShapePlugins.cpp b/DDCore/src/plugins/ShapePlugins.cpp index 0d4f629ad34c1501f34e0f270cfe986c7a8b5275..470beb7ec2a8fb49a92d4dc9ec2c3eeae9ca69c0 100644 --- a/DDCore/src/plugins/ShapePlugins.cpp +++ b/DDCore/src/plugins/ShapePlugins.cpp @@ -12,13 +12,13 @@ //========================================================================== // Framework include files -#include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/Printout.h" -#include "XML/Utilities.h" -#include "../ShapeTags.h" -#include "TGeoShapeAssembly.h" -#include "TSystem.h" -#include "TClass.h" +#include <DD4hep/DetFactoryHelper.h> +#include <DD4hep/Printout.h> +#include <XML/Utilities.h> +#include <DD4hep/ShapeTags.h> +#include <TGeoShapeAssembly.h> +#include <TSystem.h> +#include <TClass.h> // C/C++ include files #include <fstream> diff --git a/DDCore/src/plugins/StandardPlugins.cpp b/DDCore/src/plugins/StandardPlugins.cpp index 4c8d3aa4cf60ed860892523999de9af893b77779..7efef5c7fac63c1ff6e5923665d75b147d362715 100644 --- a/DDCore/src/plugins/StandardPlugins.cpp +++ b/DDCore/src/plugins/StandardPlugins.cpp @@ -13,32 +13,32 @@ // Framework include files #define DD4HEP_MUST_USE_DETECTORIMP_H -#include "DD4hep/Detector.h" -#include "DD4hep/DetectorImp.h" -#include "DD4hep/Memory.h" -#include "DD4hep/DD4hepUI.h" -#include "DD4hep/Factories.h" -#include "DD4hep/Printout.h" -#include "DD4hep/DD4hepUnits.h" -#include "DD4hep/DetectorTools.h" -#include "DD4hep/PluginCreators.h" -#include "DD4hep/VolumeProcessor.h" -#include "DD4hep/DetectorProcessor.h" -#include "DD4hep/DD4hepRootPersistency.h" -#include "XML/DocumentHandler.h" -#include "XML/XMLElements.h" -#include "XML/XMLTags.h" +#include <DD4hep/Detector.h> +#include <DD4hep/DetectorImp.h> +#include <DD4hep/Memory.h> +#include <DD4hep/DD4hepUI.h> +#include <DD4hep/Factories.h> +#include <DD4hep/Printout.h> +#include <DD4hep/DD4hepUnits.h> +#include <DD4hep/DetectorTools.h> +#include <DD4hep/PluginCreators.h> +#include <DD4hep/VolumeProcessor.h> +#include <DD4hep/DetectorProcessor.h> +#include <DD4hep/DD4hepRootPersistency.h> +#include <XML/DocumentHandler.h> +#include <XML/XMLElements.h> +#include <XML/XMLTags.h> // ROOT includes -#include "TInterpreter.h" -#include "TGeoElement.h" -#include "TGeoManager.h" -#include "TGeoVolume.h" -#include "TSystem.h" -#include "TClass.h" -#include "TRint.h" +#include <TInterpreter.h> +#include <TGeoElement.h> +#include <TGeoManager.h> +#include <TGeoVolume.h> +#include <TSystem.h> +#include <TClass.h> +#include <TRint.h> #if ROOT_VERSION_CODE >= ROOT_VERSION(6,17,0) -#include "TGDMLMatrix.h" +#include <TGDMLMatrix.h> #endif // C/C++ include files diff --git a/DDDigi/CMakeLists.txt b/DDDigi/CMakeLists.txt index bac993b915670d6881c18daaeeae044d5a007924..e6fb109605b4b7aa9cf8e88bafe338d8945637a2 100644 --- a/DDDigi/CMakeLists.txt +++ b/DDDigi/CMakeLists.txt @@ -14,15 +14,24 @@ file(GLOB DDDigi_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/src/*.cpp add_library(DDDigi ${DDDigi_SOURCES}) add_library(DD4hep::DDDigi ALIAS DDDigi) +## FIND_PACKAGE(GSL REQUIRED QUIET) +if(GSL_FOUND) + dd4hep_print( "|++> GSL_INCLUDE_DIR --> ${GSL_INCLUDE_DIR}") + dd4hep_print( "|++> GSL_LIBRARY --> ${GSL_LIBRARY}") + dd4hep_print( "|++> GSL found. DDDigi will run multi threaded.") +endif() + target_link_libraries(DDDigi PUBLIC - DD4hep::DDCore Boost::boost ROOT::Core ROOT::Geom ROOT::GenVector ROOT::RIO) + DD4hep::DDCore Boost::boost ROOT::Core ROOT::Geom ROOT::GenVector ROOT::RIO ${GSL_LIBRARY}) target_include_directories(DDDigi PUBLIC + ${GSL_INCLUDE_DIR} $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include> $<INSTALL_INTERFACE:include> ) + FIND_PACKAGE(TBB QUIET) if(TBB_FOUND) dd4hep_print( "|++> TBB_INCLUDE_DIR --> ${TBB_INCLUDE_DIR}") diff --git a/DDDigi/include/DDDigi/DigiGaussianNoise.h b/DDDigi/include/DDDigi/DigiGaussianNoise.h new file mode 100644 index 0000000000000000000000000000000000000000..8b36d807ea509c6b902e7999afa9c2a2ab1a69e6 --- /dev/null +++ b/DDDigi/include/DDDigi/DigiGaussianNoise.h @@ -0,0 +1,62 @@ +//========================================================================== +// AIDA Detector description implementation +//-------------------------------------------------------------------------- +// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN) +// All rights reserved. +// +// For the licensing terms see $DD4hepINSTALL/LICENSE. +// For the list of contributors see $DD4hepINSTALL/doc/CREDITS. +// +// Author : M.Frank +// +//========================================================================== +#ifndef DD4HEP_DDDIGI_DIGIGAUSSIANNOISE_H +#define DD4HEP_DDDIGI_DIGIGAUSSIANNOISE_H + +/// Framework include files +#include "DDDigi/DigiSignalProcessor.h" + +/// Namespace for the AIDA detector description toolkit +namespace dd4hep { + + /// Namespace for the Digitization part of the AIDA detector description toolkit + namespace digi { + + // Forward declarations + class DigiGaussianNoise; + + /// Generic noise source with a gaussian distribution + /** + * Generate gaussian noise as it appears e.g. from electronic noise + * with a given mean and a given sigma + * + * \author M.Frank + * \version 1.0 + * \ingroup DD4HEP_DIGITIZATION + */ + class DigiGaussianNoise : public DigiSignalProcessor { + protected: + /// Property: Mean value of the + double m_mean = 0.0; + /// Property: Variance of the energy distribution in electron Volt. MANDATORY! + double m_sigma = -1.0; + /// Property: Cut-off parameter. Do nothing if existing energy deposit is below threshold + double m_cutoff = -1.0; + + protected: + /// Define standard assignments and constructors + DDDIGI_DEFINE_ACTION_CONSTRUCTORS(DigiGaussianNoise); + + public: + /// Standard constructor + DigiGaussianNoise(const DigiKernel& kernel, const std::string& nam); + /// Default destructor + virtual ~DigiGaussianNoise(); + /// Initialize the noise source + virtual void initialize() override; + /// Callback to read event gaussiannoise + virtual double operator()(const DigiCellData& data) const override; + }; + } // End namespace digi +} // End namespace dd4hep +#endif // DD4HEP_DDDIGI_DIGIGAUSSIANNOISE_H diff --git a/DDDigi/include/DDDigi/DigiRandomNoise.h b/DDDigi/include/DDDigi/DigiRandomNoise.h new file mode 100644 index 0000000000000000000000000000000000000000..8527796a6277c42f5e97833eab8c1c8a03326a33 --- /dev/null +++ b/DDDigi/include/DDDigi/DigiRandomNoise.h @@ -0,0 +1,71 @@ +//========================================================================== +// AIDA Detector description implementation +//-------------------------------------------------------------------------- +// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN) +// All rights reserved. +// +// For the licensing terms see $DD4hepINSTALL/LICENSE. +// For the list of contributors see $DD4hepINSTALL/doc/CREDITS. +// +// Author : M.Frank +// +//========================================================================== +#ifndef DD4HEP_DDDIGI_DIGIRANDOMNOISE_H +#define DD4HEP_DDDIGI_DIGIRANDOMNOISE_H + +/// Framework include files +#include "DDDigi/DigiSignalProcessor.h" +#include "DDDigi/FalphaNoise.h" + +/// Namespace for the AIDA detector description toolkit +namespace dd4hep { + + /// Namespace for the Digitization part of the AIDA detector description toolkit + namespace digi { + + // Forward declarations + class DigiRandomNoise; + + /// Generic noise source for colored noise: white, pink and brown + /** + * Generic noise source for colored noise: white, pink and brown + * Uses internall a 1 / f**alpha noise generator with + * alpha = 0 -> white noise + * alpha = 1 -> pink noise (decay 10 db per decade) + * alpha = 2 -> red (brownian) noise (decay 20 db per decade) + * See https://en.wikipedia.org/wiki/White_noise for details on colored noise. + * See the header FalphaNoise.h and references therein for details + * about the generation. + * + * \author M.Frank + * \version 1.0 + * \ingroup DD4HEP_DIGITIZATION + */ + class DigiRandomNoise : public DigiSignalProcessor { + protected: + /// Property: Alpha parameter of the distribution + double m_alpha = 1.0; + /// Property: Variance of the energy distribution in electron Volt. MANDATORY! + double m_variance = -1; + /// Property: Number of IRR poles for the noise generator (5 should fit nearly everything) + double m_poles = 5; + + /// Noise generator + detail::FalphaNoise m_noise; + protected: + /// Define standard assignments and constructors + DDDIGI_DEFINE_ACTION_CONSTRUCTORS(DigiRandomNoise); + + public: + /// Standard constructor + DigiRandomNoise(const DigiKernel& kernel, const std::string& nam); + /// Default destructor + virtual ~DigiRandomNoise(); + /// Initialize the noise source + virtual void initialize() override; + /// Callback to read event randomnoise + virtual double operator()(const DigiCellData& data) const; + }; + } // End namespace digi +} // End namespace dd4hep +#endif // DD4HEP_DDDIGI_DIGIRANDOMNOISE_H diff --git a/DDDigi/include/DDDigi/DigiSignalProcessor.h b/DDDigi/include/DDDigi/DigiSignalProcessor.h index fc7298d46a5f9507fefd741d8d31e27689524cde..782d6318275837fdea81fdab5d261f0d1d385115 100644 --- a/DDDigi/include/DDDigi/DigiSignalProcessor.h +++ b/DDDigi/include/DDDigi/DigiSignalProcessor.h @@ -37,17 +37,19 @@ namespace dd4hep { */ class DigiSignalProcessor : public DigiAction { protected: + /// Flag to check if initialized was called + bool m_initialized = false; + /// Define standard assignments and constructors DDDIGI_DEFINE_ACTION_CONSTRUCTORS(DigiSignalProcessor); - /// Main functional callback - virtual void execute(DigiContext&) const {} - public: /// Standard constructor DigiSignalProcessor(const DigiKernel& kernel, const std::string& nam); /// Default destructor virtual ~DigiSignalProcessor(); + /// Initialize the noise source + virtual void initialize(); /// Callback to read event signalprocessor virtual double operator()(const DigiCellData& data) const = 0; }; diff --git a/DDDigi/include/DDDigi/FalphaNoise.h b/DDDigi/include/DDDigi/FalphaNoise.h new file mode 100644 index 0000000000000000000000000000000000000000..e317f1b8ad4769d2efc5844f8b72b050701b49c6 --- /dev/null +++ b/DDDigi/include/DDDigi/FalphaNoise.h @@ -0,0 +1,174 @@ +//========================================================================== +// AIDA Detector description implementation +//-------------------------------------------------------------------------- +// +// This code is originally written by Sampo Niskanen in java. +// Translated to C++ by M.Frank +// +// See original copyright notice below +//========================================================================== +// +// PinkNoise.java - a pink noise generator +// +// Copyright (c) 2008, Sampo Niskanen <sampo.niskanen@iki.fi> +// All rights reserved. +// Source: http://www.iki.fi/sampo.niskanen/PinkNoise/ +// +// Distrubuted under the BSD license: +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions +// are met: +// +// - Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// - Redistributions in binary form must reproduce the above +// copyright notice, this list of conditions and the following +// disclaimer in the documentation and/or other materials provided +// with the distribution. +// +// - Neither the name of the copyright owners nor contributors may be +// used to endorse or promote products derived from this software +// without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS +// FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE +// COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, +// INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, +// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +// LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +// LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN +// ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// +//========================================================================== + +#ifndef DD4HEP_DDDIGI_FALPHANOISE_H +#define DD4HEP_DDDIGI_FALPHANOISE_H + +/// C/C++ include files +#include <vector> +#include <random> + +/// Namespace for the AIDA detector description toolkit +namespace dd4hep { + + /// Namespace for the Digitization part of the AIDA detector description toolkit + namespace detail { + + /// Source of pink noise with a power spectrum proportional to 1/f^alpha. + /** + * A class that provides a source of pink noise with a power spectrum + * density (PSD) proportional to 1/f^alpha. "Regular" pink noise has a + * PSD proportional to 1/f, i.e. alpha=1. However, many natural + * systems may require a different PSD proportionality. The value of + * alpha may be from 0 to 2, inclusive. The special case alpha=0 + * results in white noise (directly generated random numbers) and + * alpha=2 results in brown noise (integrated white noise). + * <p> + * The values are computed by applying an IIR filter to generated + * Gaussian random numbers. The number of poles used in the filter + * may be specified. For each number of poles there is a limiting + * frequency below which the PSD becomes constant. Values as low as + * 1-3 poles produce relatively good results, however these values + * will be concentrated near zero. Using a larger number of poles + * will allow more low frequency components to be included, leading to + * more variation from zero. However, the sequence is stationary, + * that is, it will always return to zero even with a large number of + * poles. + * <p> + * The distribution of values is very close to Gaussian with mean + * zero, but the variance depends on the number of poles used. + * The resulting distribution is almost Gaussian, but has a relatively + * larger amount of large values. To re-normalize the variance a + * calibration method normalize() is provided. + * <p> + * The IIR filter used by this class is presented by N. Jeremy Kasdin, + * Proceedings of the IEEE, Vol. 83, No. 5, May 1995, p. 822. + * + * \author Sampo Niskanen <sampo.niskanen@iki.fi> + * \author M.Frank + * \version 1.0 + * \ingroup DD4HEP_DIGITIZATION + */ + class FalphaNoise { + public: + /// Minimal class to mimic a std::random_engine + struct random_engine_wrapper { + typedef unsigned int result_type; + static constexpr result_type min() { return std::default_random_engine::min(); } + static constexpr result_type max() { return std::default_random_engine::max(); } + static constexpr double range() { return double(max()) - double(min()); } + virtual result_type operator()() const = 0; + }; + /// Concrete wrapper of a user defined random engine + template <typename ENGINE> struct random_engine : public random_engine_wrapper { + ENGINE& engine; + random_engine(ENGINE& e) : engine(e) {} + virtual result_type operator()() const { + static constexpr double norm = random_engine_wrapper::range() / (double(ENGINE::max() - ENGINE::min())); + return result_type(double(engine()) * norm); + } + }; + + protected: +#ifdef __GSL_FALPHA_NOISE + long ptr = -1; + double * arr = nullptr; +#endif + /// Number of poles used for the approximation + size_t m_poles = 0; + /// The variance of the pink noise distribution + double m_variance = -1e0; + /// The alpha coefficient for the 1/f**alpha generated noise + double m_alpha = -1e0; + /// Array of multipliers + std::vector<double> m_multipliers {}; + /// Array of buffered values + std::vector<double> m_values {}; + /// Normal distribution with variance for the white noise generation + std::normal_distribution<double> m_distribution; + + /// Internal usage for the transformation from the normal distruted white noise to pink noise + double compute(double gaussian); + /// Internal implementation of the normalization of the variance + void normalizeVariance(const random_engine_wrapper& engine, size_t shots); + + public: + /// Default constructor. Requires a later call to init(...) + FalphaNoise() = default; + /// Initializing constructor. Leaves the enerator full functional + FalphaNoise(size_t poles, double alpha, double variance); + /// Initializing constructor with 5 poles. Leaves the enerator full functional + FalphaNoise(double alpha, double variance); + /// Default destructor + virtual ~FalphaNoise() = default; + /// Access variance value + double variance() const { return m_variance; } + /// Access alpha value + double alpha() const { return m_alpha; } + /// Initialize the 1/f**alpha random generator. If already called reconfigures. + void init(size_t poles, double alpha, double variance); + /// Approximatively compute proper variance and apply computed value + void normalize(size_t shots=10000); + /// Approximatively compute proper variance using external ranfom engine and apply computed value + template <typename ENGINE> void normalize(ENGINE& engine, size_t shots=10000); + /// Retrieve the next random number of the sequence + template <typename ENGINE> double operator()(ENGINE& engine); + }; + + /// Retrieve the next random number of the sequence + template <typename ENGINE> inline double FalphaNoise::operator()(ENGINE& engine) { + return compute(m_distribution(engine)); + } + /// Approximatively compute proper variance using external ranfom engine and apply computed value + template <typename ENGINE> void FalphaNoise::normalize(ENGINE& engine, size_t shots) { + normalizeVariance(random_engine<ENGINE>(engine), shots); + } + } // End namespace detail +} // End namespace dd4hep +#endif // DD4HEP_DDDIGI_FALPHANOISE_H diff --git a/DDDigi/src/DigiGaussianNoise.cpp b/DDDigi/src/DigiGaussianNoise.cpp new file mode 100644 index 0000000000000000000000000000000000000000..054bfd2f2f704f9a7f3cd4e29553bc8709766251 --- /dev/null +++ b/DDDigi/src/DigiGaussianNoise.cpp @@ -0,0 +1,41 @@ +//========================================================================== +// AIDA Detector description implementation +//-------------------------------------------------------------------------- +// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN) +// All rights reserved. +// +// For the licensing terms see $DD4hepINSTALL/LICENSE. +// For the list of contributors see $DD4hepINSTALL/doc/CREDITS. +// +// Author : M.Frank +// +//========================================================================== + +// Framework include files +#include "DD4hep/InstanceCount.h" +#include "DDDigi/DigiGaussianNoise.h" + +using namespace dd4hep::digi; + +/// Standard constructor +DigiGaussianNoise::DigiGaussianNoise(const DigiKernel& krnl, const std::string& nam) + : DigiSignalProcessor(krnl, nam) +{ + InstanceCount::increment(this); +} + +/// Default destructor +DigiGaussianNoise::~DigiGaussianNoise() { + InstanceCount::decrement(this); +} + +/// Initialize the noise source +void DigiGaussianNoise::initialize() { + + DigiSignalProcessor::initialize(); +} + +/// Callback to read event gaussiannoise +double DigiGaussianNoise::operator()(const DigiCellData& /* data */) const { + return 0.0; +} diff --git a/DDDigi/src/DigiRandomNoise.cpp b/DDDigi/src/DigiRandomNoise.cpp new file mode 100644 index 0000000000000000000000000000000000000000..fcfd38675d9b01dbe64d2e9f3c53df6e4b6ececf --- /dev/null +++ b/DDDigi/src/DigiRandomNoise.cpp @@ -0,0 +1,43 @@ +//========================================================================== +// AIDA Detector description implementation +//-------------------------------------------------------------------------- +// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN) +// All rights reserved. +// +// For the licensing terms see $DD4hepINSTALL/LICENSE. +// For the list of contributors see $DD4hepINSTALL/doc/CREDITS. +// +// Author : M.Frank +// +//========================================================================== + +// Framework include files +#include "DD4hep/InstanceCount.h" +#include "DDDigi/DigiRandomNoise.h" + +using namespace dd4hep::digi; + +/// Standard constructor +DigiRandomNoise::DigiRandomNoise(const DigiKernel& krnl, const std::string& nam) + : DigiSignalProcessor(krnl, nam) +{ + InstanceCount::increment(this); +} + +/// Default destructor +DigiRandomNoise::~DigiRandomNoise() { + InstanceCount::decrement(this); +} + +/// Initialize the noise source +void DigiRandomNoise::initialize() { + std::default_random_engine generator; + m_noise.init(m_poles, m_alpha, m_variance); + m_noise.normalize(generator, 5000); + DigiSignalProcessor::initialize(); +} + +/// Callback to read event randomnoise +double DigiRandomNoise::operator()(const DigiCellData& /* data */) const { + return 0.0; +} diff --git a/DDDigi/src/DigiSignalProcessor.cpp b/DDDigi/src/DigiSignalProcessor.cpp index 3595541e05ef61db395ed9f9dac178790b1e4243..d24be895c55559687b0d16eba791a9f094e7e55b 100644 --- a/DDDigi/src/DigiSignalProcessor.cpp +++ b/DDDigi/src/DigiSignalProcessor.cpp @@ -27,3 +27,8 @@ dd4hep::digi::DigiSignalProcessor::~DigiSignalProcessor() { InstanceCount::decrement(this); } +/// Initialize the noise source +void dd4hep::digi::DigiSignalProcessor::initialize() { + m_initialized = true; +} + diff --git a/DDDigi/src/DigiSubdetectorSequence.cpp b/DDDigi/src/DigiSubdetectorSequence.cpp index 7476aec8d7920128ee3c6421587e5df447199a99..e49f4d3da226ba27f40832e523d8f22144f405ba 100644 --- a/DDDigi/src/DigiSubdetectorSequence.cpp +++ b/DDDigi/src/DigiSubdetectorSequence.cpp @@ -44,7 +44,7 @@ DigiSubdetectorSequence::~DigiSubdetectorSequence() { InstanceCount::decrement(this); } -/// Iniitalize subdetector sequencer +/// Initialize subdetector sequencer void DigiSubdetectorSequence::initialize() { info("Initializing detector sequencer for detector: %s",m_detectorName.c_str()); m_detector = subdetector(m_detectorName); diff --git a/DDDigi/src/FalphaNoise.cpp b/DDDigi/src/FalphaNoise.cpp new file mode 100644 index 0000000000000000000000000000000000000000..40d085c2cac54c694e3efce12e7536343d991599 --- /dev/null +++ b/DDDigi/src/FalphaNoise.cpp @@ -0,0 +1,426 @@ +//========================================================================== +// AIDA Detector description implementation +//-------------------------------------------------------------------------- +// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN) +// All rights reserved. +// +// For the licensing terms see $DD4hepINSTALL/LICENSE. +// For the list of contributors see $DD4hepINSTALL/doc/CREDITS. +// +// Author : M.Frank +// +//========================================================================== + +/// Framework include files +#include <DDDigi/FalphaNoise.h> + +/// C/C++ include files +#include <cmath> +#include <iostream> +#include <stdexcept> + +using namespace std; +using namespace dd4hep::detail; + +#ifdef __GSL_FALPHA_NOISE +#ifndef __CNOISE_C +#define __CNOISE_C + +/* +* This code is distributed under GPLv3 +* +* This code generates correlated or colored noise with a 1/f^alpha power +* spectrum distribution. +* +* It uses the algorithm by: +* +* Jeremy Kasdin +* Discrete Simulation of Colored Noise and Stochastic Processes and $ 1/f^\alpha $ Power Law Noise Generation +* Proceedings of the IEEE +* Volume 83, Number 5, 1995, pages 802-827. +* +* This code uses GSL fast Fourier transform gsl_fft_complex_forward(...) +* and the GCC rand() functions +* +* Code Author: Miroslav Stoyanov, Jan 2011 +* +* Copyright (C) 2011 Miroslav Stoyanov +* +* 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. +* +* Since the GNU General Public License is longer than this entire code, +* a copy of it can be obtained separately at <http://www.gnu.org/licenses/> +* +* updated June 2011: fixed a small bug causing "nan" to be returned sometimes +* +*/ + +//#include "cnoise.h" +#ifndef __CNOISE_H +#define __CNOISE_H + +/* +* This code isdistrubuted under GPLv3 +* +* This code generates corelated or colored noise with 1/f^alpha power spectrum distribution. +* It uses the algorithm by: +* +* Jeremy Kasdin +* Discrete Simulation of Colored Noise and Stochastic Processes and $ 1/f^\alpha $ Power Law Noise Generation +* Proceedings of the IEEE +* Volume 83, Number 5, 1995, pages 802-827. +* +* This code uses GSL fast Fourier transform gsl_fft_complex_forward(...) and the GCC rand() functions +* +* Code Author: Miroslav Stoyanov, Jan 2011 +* +* Copyright (C) 2011 Miroslav Stoyanov +* +* 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. +* +* Since the GNU General Public License is longer than this entire code, +* a copy of it can be obtained separately at <http://www.gnu.org/licenses/> +* +*/ + +#include <stdlib.h> +#include <stdio.h> +#include <math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_fft_complex.h> + +double cnoise_uniform01 ( ); + +void cnoise_two_gaussian_truncated ( double *a, double *b, double std, double range ); + +void cnoise_generate_colored_noise ( double *x, int N, double alpha, double std ); + +// generates a vector x of size N with a 1/f^alpha frequency distribution +// std is the standard deviation of the underlying gaussian distribution +// Variables: double *x - Input: allocated vector of size N +// Output: a realization of 1/f^alpha noise vector of size N, using an undelying Gaussian (0,std) +// int N - Input: the size of the vector +// double alpha - Input: the decay rate for the power spectrum (i.e. 1/f^alpha +// double std - Input: the standard deviation of the underlying Gaussian distribution +// NOTE: you should call srand( seed ) before you call dcnoise_generate_colored_noise + +void cnoise_generate_colored_noise_uniform( double *x, int N, double alpha, double range ); +// same as cnoise_generate_colored_noise_uniform, except the white noise vector comes from +// uniform distribution on (-range,+range) + +void cnoise_generate_colored_noise_truncated( double *x, int N, double alpha, double std, double range ); +// same as cnoise_generate_colored_noise_uniform, except the white noise vector comes from +// truncated Gaussian distribution with mean 0, standard deviation std and truncated to (-range,range) + + + +#endif + + +#define _CNOISE_PI 3.14159265358979323846 // pi + +/******************************************************************************/ + +double cnoise_uniform01 ( ) + +/******************************************************************************/ +{ + return ( ( (double) rand() + 1.0 ) / ( (double) RAND_MAX + 1.0 ) ); +}; +/******************************************************************************/ + +void cnoise_two_gaussian_truncated ( double *a, double *b, double std, double range ) + +/******************************************************************************/ +{ + double gauss_u = cnoise_uniform01(); + double gauss_v = cnoise_uniform01(); + *a = std * sqrt( - 2 * log( gauss_u ) ) * cos( 2 * _CNOISE_PI * gauss_v ); + *b = std * sqrt( - 2 * log( gauss_u ) ) * sin( 2 * _CNOISE_PI * gauss_v ); + while ( (*a < -range) || (*a > range) ){ + gauss_u = cnoise_uniform01(); gauss_v = cnoise_uniform01(); + *a = std * sqrt( - 2 * log( gauss_u ) ) * cos( 2 * _CNOISE_PI * gauss_v ); + }; + while ( (*b < -range) || (*b > range) ){ + gauss_u = cnoise_uniform01(); gauss_v = cnoise_uniform01(); + *b = std * sqrt( - 2 * log( gauss_u ) ) * cos( 2 * _CNOISE_PI * gauss_v ); + }; +}; +/******************************************************************************/ + +void cnoise_generate_colored_noise ( double *x, int N, double alpha, double std ) + +/******************************************************************************/ +{ + int i; + double tmp; + double gauss_u, gauss_v; + double * fh = (double*)malloc( 4*N*sizeof(double) ); + double * fw = (double*)malloc( 4*N*sizeof(double) ); + + fh[0] = 1.0; fh[1] = 0.0; + for( i=1; i<N; i++ ){ + fh[2*i] = fh[2*(i-1)] * ( 0.5 * alpha + (double)(i-1) ) / ((double) i ); + fh[2*i+1] = 0.0; + }; + for( i=0; i<N; i+=2 ){ + gauss_u = cnoise_uniform01(); + gauss_v = cnoise_uniform01(); + fw[2*i] = std * sqrt( - 2 * log( gauss_u ) ) * cos( 2 * _CNOISE_PI * gauss_v ); + fw[2*i+1] = 0.0; + fw[2*i+2] = std * sqrt( - 2 * log( gauss_u ) ) * sin( 2 * _CNOISE_PI * gauss_v ); + fw[2*i+3] = 0.0; + }; + for( i=2*N; i<4*N; i++ ){ + fh[i] = fw[i] = 0.0; + }; + + gsl_fft_complex_wavetable* wavetable = gsl_fft_complex_wavetable_alloc(2*N); + gsl_fft_complex_workspace* workspace = gsl_fft_complex_workspace_alloc(2*N); + + gsl_fft_complex_forward (fh, 1, 2*N, wavetable, workspace); + gsl_fft_complex_forward (fw, 1, 2*N, wavetable, workspace); + + for( i=0; i<N+1; i++ ){ + tmp = fw[2*i]; + fw[2*i] = tmp*fh[2*i] - fw[2*i+1]*fh[2*i+1]; + fw[2*i+1] = tmp*fh[2*i+1] + fw[2*i+1]*fh[2*i]; + }; + + fw[0] /= 2.0; fw[1] /= 2.0; + fw[2*N] /= 2.0; fw[2*N+1] /= 2.0; + + for( i=N+1; i<2*N; i++ ){ + fw[2*i] = 0.0; fw[2*i+1] = 0.0; + }; + + gsl_fft_complex_inverse( fw, 1, 2*N, wavetable, workspace); + + for( i=0; i<N; i++ ){ + x[i] = 2.0*fw[2*i]; + }; + + free( fh ); + free( fw ); + + gsl_fft_complex_wavetable_free (wavetable); + gsl_fft_complex_workspace_free (workspace); +}; +/******************************************************************************/ + +void cnoise_generate_colored_noise_uniform( double *x, int N, double alpha, + double range ) + +/******************************************************************************/ +{ + gsl_fft_complex_wavetable * wavetable; + gsl_fft_complex_workspace * workspace; + + int i; + double tmp; + double * fh = (double*)malloc( 4*N*sizeof(double) ); + double * fw = (double*)malloc( 4*N*sizeof(double) ); + + fh[0] = 1.0; fh[1] = 0.0; + fw[0] = 2.0*range*cnoise_uniform01()-range; fw[1] = 0.0; + for( i=1; i<N; i++ ){ + fh[2*i] = fh[2*(i-1)] * ( 0.5 * alpha + (double)(i-1) ) / ((double) i ); + fh[2*i+1] = 0.0; + fw[2*i] = 2.0*range*cnoise_uniform01()-range; fw[2*i+1] = 0.0; + }; + for( i=2*N; i<4*N; i++ ){ fh[i] = 0.0; fw[i] = 0.0; }; + + wavetable = gsl_fft_complex_wavetable_alloc(2*N); + workspace = gsl_fft_complex_workspace_alloc(2*N); + + gsl_fft_complex_forward (fh, 1, 2*N, wavetable, workspace); + gsl_fft_complex_forward (fw, 1, 2*N, wavetable, workspace); + + for( i=0; i<N+1; i++ ){ + tmp = fw[2*i]; + fw[2*i] = tmp*fh[2*i] - fw[2*i+1]*fh[2*i+1]; + fw[2*i+1] = tmp*fh[2*i+1] + fw[2*i+1]*fh[2*i]; + }; + + fw[0] /= 2.0; fw[1] /= 2.0; + fw[2*N] /= 2.0; fw[2*N+1] /= 2.0; + + for( i=N+1; i<2*N; i++ ){ + fw[2*i] = 0.0; fw[2*i+1] = 0.0; + }; + + gsl_fft_complex_inverse( fw, 1, 2*N, wavetable, workspace); + + for( i=0; i<N; i++ ){ + x[i] = 2.0*fw[2*i]; + }; + + free( fh ); + free( fw ); + + gsl_fft_complex_wavetable_free (wavetable); + gsl_fft_complex_workspace_free (workspace); +}; +/******************************************************************************/ + +void cnoise_generate_colored_noise_truncated( double *x, int N, double alpha, double std, double range ) + +/******************************************************************************/ +{ + gsl_fft_complex_wavetable * wavetable; + gsl_fft_complex_workspace * workspace; + + int i; + double tmp; + double * fh = (double*)malloc( 4*N*sizeof(double) ); + double * fw = (double*)malloc( 4*N*sizeof(double) ); + + fh[0] = 1.0; fh[1] = 0.0; + for( i=1; i<N; i++ ){ + fh[2*i] = fh[2*(i-1)] * ( 0.5 * alpha + (double)(i-1) ) / ((double) i ); + fh[2*i+1] = 0.0; + }; + for( i=0; i<N; i+=2 ){ + cnoise_two_gaussian_truncated( &fw[2*i], &fw[2*i+2], std, range ); + }; + for( i=2*N; i<4*N; i++ ){ fh[i] = 0.0; fw[i] = 0.0; }; + + wavetable = gsl_fft_complex_wavetable_alloc(2*N); + workspace = gsl_fft_complex_workspace_alloc(2*N); + + gsl_fft_complex_forward (fh, 1, 2*N, wavetable, workspace); + gsl_fft_complex_forward (fw, 1, 2*N, wavetable, workspace); + + for( i=0; i<N+1; i++ ){ + tmp = fw[2*i]; + fw[2*i] = tmp*fh[2*i] - fw[2*i+1]*fh[2*i+1]; + fw[2*i+1] = tmp*fh[2*i+1] + fw[2*i+1]*fh[2*i]; + }; + + fw[0] /= 2.0; fw[1] /= 2.0; + fw[2*N] /= 2.0; fw[2*N+1] /= 2.0; + + for( i=N+1; i<2*N; i++ ){ + fw[2*i] = 0.0; fw[2*i+1] = 0.0; + }; + + gsl_fft_complex_inverse( fw, 1, 2*N, wavetable, workspace); + + for( i=0; i<N; i++ ){ + x[i] = 2.0*fw[2*i]; + }; + + free ( fh ); + free ( fw ); + + gsl_fft_complex_wavetable_free (wavetable); + gsl_fft_complex_workspace_free (workspace); +}; + +#endif +#endif + +/// Initializing constructor. Leaves the enerator full functional +FalphaNoise::FalphaNoise(size_t poles, double alpha, double variance) { + init(poles, alpha, variance); +} + +/// Initializing constructor with 5 poles. Leaves the enerator full functional +FalphaNoise::FalphaNoise(double alpha, double variance) { + init(5, alpha, variance); +} + +/// Initialize the 1/f**alpha random generator. If already called reconfigures. +void FalphaNoise::init(size_t poles, double alpha, double variance) { + m_poles = poles; + m_alpha = alpha; + m_variance = variance; + if ( alpha < 0e0 || alpha > 2e0 ) { + throw std::runtime_error("FalphaNoise: Invalid value for alpha: must be inside the interval [0,2]"); + } + else if ( variance < 1e-10 ) { + throw std::runtime_error("FalphaNoise: Invalid variance. Must be bigger than NULL!"); + } + m_values.clear(); + m_multipliers.clear(); + m_distribution = std::normal_distribution<double>(0.0, m_variance); + m_values.resize(m_poles+1,0e0); + m_multipliers.reserve(m_poles); + + double val = 1e0; + for ( size_t i=0; i<m_poles; ++i) { + val *= (double(i) - alpha/2e0) / double(i+1); + m_multipliers.emplace_back(val); + } + // Fill the history with random values + std::default_random_engine generator; + for ( size_t i=0; i < 5*m_poles; i++) + (*this)(generator); +} + +/// Approximative compute proper variance and apply computed value +void FalphaNoise::normalize(size_t shots) { + std::default_random_engine engine; + normalize(engine, shots); +} + +/// Internal usage to normalize using user defined random engine +void FalphaNoise::normalizeVariance(const random_engine_wrapper& generator, size_t shots) { + double mean = 0e0, mean2 = 0e0, var; + auto tmp = m_multipliers; + double val = 1e0; + /// We have to correct for the case of white noise: alpha=0 + for ( size_t i=0; i<m_poles; ++i) { + val *= (double(i) - 0.5) / double(i+1); + m_multipliers[i] = val; + } + for ( size_t i=0; i < shots; i++) { + var = (*this)(generator); + mean += var; + mean2 += var*var; + } + mean /= double(shots); + var = std::sqrt(mean2/double(shots) - mean*mean); + m_variance *= m_variance/var; + m_multipliers = tmp; + m_distribution = std::normal_distribution<double>(0.0, m_variance); +} + +/// Retrieve the next random number of the sequence +double FalphaNoise::compute(double rndm_value) { +#ifdef __GSL_FALPHA_NOISE + if ( arr == nullptr ) { + arr = new double[1000]; + } + if ( ptr < 0 ) { + cout << "Alpha: " << m_alpha << " sigma: " << m_variance << endl; + cnoise_generate_colored_noise (arr, 1000, m_alpha, m_variance); + ptr = 999; + } + --ptr; + return arr[ptr+1]; +#else + for ( int i=m_poles-1; i >= 0; --i ) { + rndm_value -= m_multipliers[i] * m_values[i]; + m_values[i+1] = m_values[i]; // highest index always drops off... + } + m_values[0] = rndm_value; + return rndm_value; +#endif +} diff --git a/examples/DDDigi/CMakeLists.txt b/examples/DDDigi/CMakeLists.txt index e6b7cf089456082b2919e95bec4099d93c50f76b..d660e26c1f705344976f65d1e8cd93baa2dacbdd 100644 --- a/examples/DDDigi/CMakeLists.txt +++ b/examples/DDDigi/CMakeLists.txt @@ -24,13 +24,13 @@ dd4hep_configure_output () set(DDDigiexamples_INSTALL ${CMAKE_INSTALL_PREFIX}/examples/DDDigi) # dd4hep_add_plugin(DDDigiExampleLib SOURCES src/*.cpp ) -target_link_libraries(DDDigiExampleLib DD4hep::DDDigi Boost::boost ROOT::Geom ROOT::GenVector ROOT::RIO) +target_link_libraries(DDDigiExampleLib DD4hep::DDDigi Boost::boost ROOT::Geom ROOT::GenVector ROOT::RIO ROOT::Gui ROOT::Hist) install(TARGETS DDDigiExampleLib LIBRARY DESTINATION lib) install(DIRECTORY scripts DESTINATION ${DDDigiexamples_INSTALL}) # dd4hep_configure_scripts (DDDigi DEFAULT_SETUP WITH_TESTS) # -# Test HepMC input reader +# Test basic processing chain dd4hep_add_test_reg(DDDigi_framework COMMAND "${CMAKE_INSTALL_PREFIX}/bin/run_test_DDDigi.sh" EXEC_ARGS python ${DDDigiexamples_INSTALL}/scripts/TestFramework.py @@ -38,4 +38,12 @@ dd4hep_add_test_reg(DDDigi_framework REGEX_FAIL "Error;ERROR;Exception" ) # +# Test colored noise factory +dd4hep_add_test_reg(DDDigi_colored_noise + COMMAND "${CMAKE_INSTALL_PREFIX}/bin/run_test_DDDigi.sh" + EXEC_ARGS geoPluginRun -ui -plugin DD4hep_FalphaNoise -shots 1000000 -variance 1 -alpha 0.5 + REGEX_PASS "FalphaNoise INFO Distribution RMS 1.0" + REGEX_FAIL "Error;ERROR;Exception" + ) +# diff --git a/examples/DDDigi/src/DigiTestFalphaNoise.cpp b/examples/DDDigi/src/DigiTestFalphaNoise.cpp new file mode 100644 index 0000000000000000000000000000000000000000..cbe437742449c0ceb3a75b6a6766d07209868511 --- /dev/null +++ b/examples/DDDigi/src/DigiTestFalphaNoise.cpp @@ -0,0 +1,164 @@ +//========================================================================== +// AIDA Detector description implementation +//-------------------------------------------------------------------------- +// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN) +// All rights reserved. +// +// For the licensing terms see $DD4hepINSTALL/LICENSE. +// For the list of contributors see $DD4hepINSTALL/doc/CREDITS. +// +// Author : M.Frank +// +//========================================================================== + +/// Framework include files +#include <DD4hep/Factories.h> +#include <DD4hep/Printout.h> +#include <DDDigi/FalphaNoise.h> + +/// C/C++ include files +#include <random> +#include <iostream> + +#include <TH1.h> +#include <TCanvas.h> +#include <TPad.h> + +using namespace dd4hep; + +/// Plugin to test the pink noise generator +/** + * Factory: DD4hep_DummyPlugin + * + * \author M.Frank + * \version 1.0 + * \date 27/11/2019 + */ +static long test_FalphaNoise(Detector& , int argc, char** argv) { + using namespace dd4hep::detail; + bool draw = false; + size_t poles = 5; + size_t shots = 10; + double alpha = 0.5, variance = 1.0; + for(int i = 0; i < argc && argv[i]; ++i) { + if ( 0 == ::strncmp("-shots",argv[i],3) ) + shots = ::atol(argv[++i]); + else if ( 0 == ::strncmp("-alpha",argv[i],3) ) + alpha = ::atof(argv[++i]); + else if ( 0 == ::strncmp("-variance",argv[i],3) ) + variance = ::atof(argv[++i]); + else if ( 0 == ::strncmp("-poles",argv[i],3) ) + poles = ::atol(argv[++i]); + else if ( 0 == ::strncmp("-draw",argv[i],3) ) + draw = true; + else { + std::cout << + "Usage: -plugin DD4hep_FalphaNoise -arg [-arg] \n" + " -shots <value> Number of samples to be generated [default: 10] \n" + " -alpha <value> Parameter for the 1/f**alpha distribution \n" + " default: 0.5 \n" + " -variance <value> Distribution variance [default: 1] \n" + " -poles <value> Number of IRR poles to fill the distribution \n" + " -draw Fill and draw hoistogram with distribution \n" + "\tArguments given: " << arguments(argc,argv) << std::endl << std::flush; + ::exit(EINVAL); + } + } + + std::default_random_engine generator; + FalphaNoise noise(poles, alpha, variance); + FalphaNoise noise0(poles, 0, variance); + FalphaNoise noise1(poles, 1, variance); + FalphaNoise noise2(poles, 1.98, variance); + std::stringstream cpara; + cpara << " distribution alpha=" << alpha << " sigma=" << variance; + TH1D* hist11 = new TH1D("D11", ("1/f**alpha"+cpara.str()).c_str(), 50, -5e0*variance, 5e0*variance); + TH1D* hist12 = new TH1D("D12", ("1/f**alpha"+cpara.str()).c_str(), 50, -5e0*variance, 5e0*variance); + TH1D* hist13 = new TH1D("D13", ("1/f**alpha"+cpara.str()).c_str(), 50, -5e0*variance, 5e0*variance); + TH1D* hist14 = new TH1D("D14", ("1/f**alpha"+cpara.str()).c_str(), 50, -5e0*variance, 5e0*variance); + TH1D* hist21 = new TH1D("D21", ("log10(1/f**alpha)"+cpara.str()).c_str(), 300, 1e0, 1e2*variance); + TH1D* hist22 = new TH1D("D22", "log10(1/f**alpha) alpha=0", 300, 1e0, 1e2*variance); + TH1D* hist23 = new TH1D("D23", "log10(1/f**alpha) alpha=1", 300, 1e0, 1e2*variance); + TH1D* hist24 = new TH1D("D24", "log10(1/f**alpha) alpha=2", 300, 1e0, 1e2*variance); + + printout(INFO, "FalphaNoise", "Executing %ld shots .... variance=%.3f alpha=%.3f", + shots, noise.variance(), noise.alpha()); + for(size_t i=0; i < shots; ++i) { + double val = noise(generator); + hist11->Fill(val, 1.); + hist21->Fill(std::exp(std::abs(val)), 1.); + + val = noise0(generator); + hist12->Fill(val, 1.); + hist22->Fill(std::exp(std::abs(val)), 1.); + + val = noise1(generator); + hist13->Fill(val, 1.); + hist23->Fill(std::exp(std::abs(val)), 1.); + + val = noise2(generator); + hist14->Fill(val, 1.); + hist24->Fill(std::exp(std::abs(val)), 1.); + } + + printout(INFO, "FalphaNoise", "Distribution Mean %10.5f", hist11->GetMean()); + printout(INFO, "FalphaNoise", "Distribution Mean Uncertainty %10.5f", hist11->GetMeanError()); + printout(INFO, "FalphaNoise", "Distribution RMS %10.5f", hist11->GetRMS()); + printout(INFO, "FalphaNoise", "Distribution RMS Uncertainty %10.5f", hist11->GetRMSError()); + hist11->GetXaxis()->SetTitle("Energy [arb.units]"); + hist11->GetYaxis()->SetTitle("Counts"); + hist21->GetXaxis()->SetTitle("exp(Energy) [arb.units]"); + hist21->GetYaxis()->SetTitle("Counts"); + + hist21->Scale(1.0 / hist21->GetBinContent(1)); + hist22->Scale(1.0 / hist22->GetBinContent(1)); + hist23->Scale(1.0 / hist23->GetBinContent(1)); + hist24->Scale(1.0 / hist24->GetBinContent(1)); + if ( draw ) { + TCanvas *c = new TCanvas("c1","multigraph",700,500); + c->Divide(1,2); + c->cd(1); + + gPad->SetGrid(); + hist12->SetLineColor(kBlue); + hist12->Draw("L"); + + hist11->SetLineWidth(1); + hist11->SetLineColor(kRed); + hist11->Draw("SAME"); + + hist11->SetLineWidth(1); + hist11->SetLineColor(kRed); + hist11->Draw("LSAME"); + + hist13->SetLineColor(kMagenta); + hist13->Draw("LSAME"); + + hist14->SetLineColor(kGreen); + hist14->Draw("LSAME"); + + c->cd(2); + gPad->SetGrid(); + gPad->SetLogx(); + gPad->SetLogy(); + hist21->Draw(""); + hist21->SetLineWidth(1); + hist21->SetLineColor(kRed); + hist21->Draw("LPSAME"); + + hist24->SetLineWidth(1); + hist24->SetLineColor(kGreen); + hist24->Draw("LPSAME"); + + hist22->SetLineWidth(1); + hist22->SetLineColor(kBlue); + hist22->Draw("LPSAME"); + hist23->SetLineWidth(1); + hist23->SetLineColor(kMagenta); + hist23->Draw("LPSAME"); + gPad->Update(); + gPad->Modified(); + } + return 1; +} +DECLARE_APPLY(DD4hep_FalphaNoise,test_FalphaNoise)