From c2dd1059dd30e9b525556a583776e82f15a9281b Mon Sep 17 00:00:00 2001 From: Gordon McCann Date: Sat, 8 Oct 2022 15:21:56 -0400 Subject: [PATCH] Finally implement randomization for conversion between integers to floating points. Eliminates binning artifacts in histograms. --- SpecProject/src/SPSAnalysisStage.cpp | 15 +++-- Specter/src/CMakeLists.txt | 1 + Specter/src/Specter.h | 1 + Specter/src/Specter/Utils/RandomGenerator.h | 64 +++++++++++++++++++++ 4 files changed, 76 insertions(+), 5 deletions(-) create mode 100644 Specter/src/Specter/Utils/RandomGenerator.h diff --git a/SpecProject/src/SPSAnalysisStage.cpp b/SpecProject/src/SPSAnalysisStage.cpp index 5411d60..31c6528 100644 --- a/SpecProject/src/SPSAnalysisStage.cpp +++ b/SpecProject/src/SPSAnalysisStage.cpp @@ -84,16 +84,20 @@ namespace Specter { { sabreFlag = true; if (hit.longEnergy > iter->second.GetValue()) - iter->second.SetValue(hit.longEnergy); + iter->second.SetValue(hit.longEnergy + RandomGenerator::GetUniformFraction()); continue; } + //Note randomization of values. In general, we want to plot/calculate as floating point. + //Especially for histograms, care has to be taken to randomize the decimal to account for binning. + //For timestamps from PSD, since we reduce accuracy by 3 orders of mag, no randomization needed. + //For timestamps from PHA, need randomization on interval from 0 to digitizer sampling period. switch (hit.id) { case s_scintLeftID: { if (hit.longEnergy > scintLeft.GetValue()) - scintLeft.SetValue(hit.longEnergy); + scintLeft.SetValue(hit.longEnergy + RandomGenerator::GetUniformFraction()); break; } case s_beamIntID: @@ -104,7 +108,7 @@ namespace Specter { case s_cathodeID: { if (hit.longEnergy > cathode.GetValue()) - cathode.SetValue(hit.longEnergy); + cathode.SetValue(hit.longEnergy + RandomGenerator::GetUniformFraction()); break; } case s_delayFrontLeftID: @@ -134,13 +138,13 @@ namespace Specter { case s_anodeFrontID: { if (hit.longEnergy > anodeFront.GetValue()) - anodeFront.SetValue(hit.longEnergy); + anodeFront.SetValue(hit.longEnergy + RandomGenerator::GetUniformFraction()); break; } case s_anodeBackID: { if (hit.longEnergy > anodeBack.GetValue()) - anodeBack.SetValue(hit.longEnergy); + anodeBack.SetValue(hit.longEnergy + RandomGenerator::GetUniformFraction()); break; } } @@ -148,6 +152,7 @@ namespace Specter { //If you want to use parameters to calculate another parameter, you //need to check that the parameter is valid (set in this event)! + //Note subsequent calculations no longer need randomization. if(delayFLTime.IsValid() && delayFRTime.IsValid()) x1.SetValue((delayFLTime.GetValue() - delayFRTime.GetValue())*0.5*0.4762); diff --git a/Specter/src/CMakeLists.txt b/Specter/src/CMakeLists.txt index 5f939f5..74c4f09 100644 --- a/Specter/src/CMakeLists.txt +++ b/Specter/src/CMakeLists.txt @@ -105,6 +105,7 @@ target_sources(Specter PRIVATE Specter/Physics/Daqromancy/DYOnlineSource.cpp Specter/Utils/Functions.h Specter/Utils/Functions.cpp + Specter/Utils/RandomGenerator.h ) #ImPlot sources diff --git a/Specter/src/Specter.h b/Specter/src/Specter.h index c8f59c8..18387e8 100644 --- a/Specter/src/Specter.h +++ b/Specter/src/Specter.h @@ -35,5 +35,6 @@ #include "Specter/Utils/TestServerLayer.h" #include "Specter/Utils/Instrumentor.h" #include "Specter/Utils/Functions.h" +#include "Specter/Utils/RandomGenerator.h" #endif diff --git a/Specter/src/Specter/Utils/RandomGenerator.h b/Specter/src/Specter/Utils/RandomGenerator.h new file mode 100644 index 0000000..7a03630 --- /dev/null +++ b/Specter/src/Specter/Utils/RandomGenerator.h @@ -0,0 +1,64 @@ +/* + RandomGenerator.h + Thread-safe implementation of stl random number generator library. A single Mersenne-Twister generator is created for + each thread. The paramters for the distributions are passed in to the function call, and a random number is calculated + for the requested distribution. Note that distributions are templated so that one can specify the type of distributions, + but most distributions only support a limited range of types. The valid types are commented for each case. + + Oct 2022 -- GWM +*/ + +#ifndef RANDOM_GENERATOR_H +#define RANDOM_GENERATOR_H + +#include + +namespace Specter { + + class RandomGenerator + { + public: + + //Valid for float, double, long double + template + static T GetUniformReal(T min, T max) + { + std::uniform_real_distribution distribution(min, max); + return distribution(GetGenerator()); + } + + //Valid for any integer type (signed or unsigned) + template + static T GetUniformInt(T min, T max) + { + std::uniform_int_distribution distribution(min, max); + return distribution(GetGenerator()); + } + + //Valid for float, double, long double + template + static T GetNormal(T mean, T sigma) + { + std::normal_distribution distribution(mean, sigma); + return distribution(GetGenerator()); + } + + //This is the most common use case, so we eliminate recreation of distribution, templating. + //For randomization of decimals in conversion from integer -> floating point for histograming + static double GetUniformFraction() + { + static std::uniform_real_distribution distribution(0.0, 1.0); + return distribution(GetGenerator()); + } + + private: + static std::mt19937& GetGenerator() + { + static thread_local std::mt19937 generator((std::random_device())()); + return generator; + } + }; +} + + +#endif \ No newline at end of file