From 85133c32de512b8a36d0f05632df9dbf5d7cb297 Mon Sep 17 00:00:00 2001 From: Nicholas Carrara Date: Thu, 5 May 2022 16:16:20 -0700 Subject: [PATCH] Fixed compiler issues with legacy code (#148) * some changes * some changes * some changes * some changes * some changes * some changes * some changes * some changes * some changes * some changes * some changes * some changes * some changes * some changes * removed plots * some changes * some changes, tried to fix code warnings * added the alpha model for LAr * added the alpha model for LAr * added the alpha model for LAr * some changes * some changes * some changes * some changes * some changes * some changes * benchmark plots finished, waiting on Mike Mooney and Justin Mueller to OK the LArNEST code for integration into the main NEST branch. Next steps are to integrate the current model into LArSoft, and then implement a dE/dx model for LAr * removed .pyc files Co-authored-by: Nicholas Carrara Co-authored-by: Nicholas Carrara --- .gitignore | 5 + CMakeLists.txt | 6 + examples/CMakeLists.txt | 2 + examples/LArNEST/CMakeLists.txt | 8 + examples/LArNEST/LArNESTBenchmarks.cpp | 113 ++++ examples/LArNEST/LegacyLArNESTBenchmarks.cpp | 113 ++++ include/Detectors/LArDetector.hh | 24 + include/NEST/LArNEST.hh | 530 ++++++++++++++++ include/NEST/TestLArSpectra.hh | 54 ++ src/LArDetector.cpp | 87 +++ src/LArNEST.cpp | 612 +++++++++++++++++++ src/NEST.cpp | 6 +- src/TestLArSpectra.cpp | 10 + utils/larnest/lar_dataset.py | 255 ++++++++ utils/larnest/parameters.py | 166 +++++ utils/larnest/plot_benchmarks.py | 40 ++ utils/larnest/utils.py | 174 ++++++ 17 files changed, 2203 insertions(+), 2 deletions(-) create mode 100644 examples/LArNEST/CMakeLists.txt create mode 100644 examples/LArNEST/LArNESTBenchmarks.cpp create mode 100644 examples/LArNEST/LegacyLArNESTBenchmarks.cpp create mode 100644 include/Detectors/LArDetector.hh create mode 100644 include/NEST/LArNEST.hh create mode 100644 include/NEST/TestLArSpectra.hh create mode 100644 src/LArDetector.cpp create mode 100644 src/LArNEST.cpp create mode 100644 src/TestLArSpectra.cpp create mode 100644 utils/larnest/lar_dataset.py create mode 100644 utils/larnest/parameters.py create mode 100644 utils/larnest/plot_benchmarks.py create mode 100644 utils/larnest/utils.py diff --git a/.gitignore b/.gitignore index 74d4af84..608f2f47 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,8 @@ /fastNEST.exe nbproject/* **/.DS_Store +build/ +*.png +utils/larnest/data/ +*.pyc +*__pycache__/ diff --git a/CMakeLists.txt b/CMakeLists.txt index 637399aa..0ceca53a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -51,14 +51,20 @@ set(NEST_CORE_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/src/VDetector.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/GammaHandler.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/ValidityTests.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/LArNEST.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/LArDetector.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/TestLArSpectra.cpp ) set(NEST_CORE_HEADERS ${CMAKE_CURRENT_SOURCE_DIR}/include/NEST/NEST.hh ${CMAKE_CURRENT_SOURCE_DIR}/include/NEST/RandomGen.hh ${CMAKE_CURRENT_SOURCE_DIR}/include/NEST/TestSpectra.hh ${CMAKE_CURRENT_SOURCE_DIR}/include/Detectors/VDetector.hh + ${CMAKE_CURRENT_SOURCE_DIR}/include/Detectors/LArDetector.hh ${CMAKE_CURRENT_SOURCE_DIR}/include/NEST/GammaHandler.hh ${CMAKE_CURRENT_SOURCE_DIR}/include/NEST/ValidityTests.hh + ${CMAKE_CURRENT_SOURCE_DIR}/include/NEST/LArNEST.hh + ${CMAKE_CURRENT_SOURCE_DIR}/include/NEST/TestLArSpectra.hh ) add_library(NESTCore ${NEST_CORE_SOURCES} ${NEST_CORE_HEADERS}) diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index b133ae64..d8ce1729 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -9,6 +9,8 @@ target_link_libraries(bareNEST PUBLIC NEST::Core) set(EXAMPLES bareNEST) +add_subdirectory("LArNEST/") + # ROOT example if(BUILD_ROOT) # https://cliutils.gitlab.io/modern-cmake/examples/root-simple/ diff --git a/examples/LArNEST/CMakeLists.txt b/examples/LArNEST/CMakeLists.txt new file mode 100644 index 00000000..2db85ae6 --- /dev/null +++ b/examples/LArNEST/CMakeLists.txt @@ -0,0 +1,8 @@ +# add_executable(bareLArNEST bareLArNEST.cpp) +# target_link_libraries(bareLArNEST PUBLIC NEST::Core) + +add_executable(LegacyLArNESTBenchmarks LegacyLArNESTBenchmarks.cpp) +target_link_libraries(LegacyLArNESTBenchmarks PUBLIC NEST::Core) + +add_executable(LArNESTBenchmarks LArNESTBenchmarks.cpp) +target_link_libraries(LArNESTBenchmarks PUBLIC NEST::Core) \ No newline at end of file diff --git a/examples/LArNEST/LArNESTBenchmarks.cpp b/examples/LArNEST/LArNESTBenchmarks.cpp new file mode 100644 index 00000000..1d26e562 --- /dev/null +++ b/examples/LArNEST/LArNESTBenchmarks.cpp @@ -0,0 +1,113 @@ +/** + * @file legacyLArNEST.cpp + * @author NEST Collaboration + * @author Nicholas Carrara [nmcarrara@ucdavis.edu] + * @brief Benchmarks for LAr ER model. This program generates ... + * @version + * @date 2022-04-14 + */ +#include +#include +#include +#include +#include + +#include "LArDetector.hh" +#include "LArNEST.hh" + +int main(int argc, char* argv[]) +{ + // this sample program takes can take two arguments, + // (optional) 1) density - the density of the LAr to use + // (optional) 2) seed - a seed for the random number generator + + double density = 1.393; + uint64_t seed = 0; + if (argc > 1) { + density = atof(argv[1]); + } + if (argc > 2) { + seed = atoi(argv[2]); + } + + // set up energy steps + int num_energy_steps = 50000; + std::vector energy_vals; + double start_val = .1; + double end_val = 1000; + double step_size = (end_val - start_val)/num_energy_steps; + + for (size_t i = 0; i < num_energy_steps; i++) + { + energy_vals.emplace_back(start_val + step_size * i); + } + std::vector particle_types = { + NEST::LArInteraction::NR, NEST::LArInteraction::ER, NEST::LArInteraction::Alpha + }; + std::vector particle_type = { + "NR", "ER", "Alpha" + }; + + LArDetector* detector = new LArDetector(); + // Construct NEST class using detector object + NEST::LArNEST larnest(detector); + // Set random seed of RandomGen class + RandomGen::rndm()->SetSeed(seed); + + // Construct NEST objects for storing calculation results + NEST::LArNESTResult result; + std::ofstream output_file; + + output_file.open("benchmarks.csv"); + output_file << "type,energy,efield,TotalYield,QuantaYield,LightYield,Nph,Ne,Nex,Nion\n"; + + // iterate over electric field values + for (size_t k = 0; k < particle_types.size(); k++) + { + std::vector electric_field; + if (particle_types[k] == NEST::LArInteraction::NR) + { + electric_field.insert(electric_field.end(), { + 1, 50, 100, 200, 250, 500, 1000, 1500, 1750, 2000 + }); + } + else if (particle_types[k] == NEST::LArInteraction::ER) + { + electric_field.insert(electric_field.end(), { + 1, 100, 200, 600, 1000, 1500, 2500, 6000, 9000, 9500 + }); + } + else + { + electric_field.insert(electric_field.end(), { + 1, 50, 100, 500, 1000, 5000, 10000, 20000 + }); + } + for (size_t v = 0; v < electric_field.size(); v++) + { + // iterate over energy values + for (size_t i = 0; i < num_energy_steps; i++) + { + result = larnest.FullCalculation( + particle_types[k], + energy_vals[i], + electric_field[v], + density, + false + ); + output_file << particle_type[k] << ","; + output_file << energy_vals[i] << ","; + output_file << electric_field[v] << ","; + output_file << result.yields.TotalYield << ","; + output_file << result.yields.QuantaYield << ","; + output_file << result.yields.LightYield << ","; + output_file << result.yields.Nph << ","; + output_file << result.yields.Ne << ","; + output_file << result.yields.Nex << ","; + output_file << result.yields.Nion << "\n"; + } + } + } + output_file.close(); + return 0; +} \ No newline at end of file diff --git a/examples/LArNEST/LegacyLArNESTBenchmarks.cpp b/examples/LArNEST/LegacyLArNESTBenchmarks.cpp new file mode 100644 index 00000000..cbc8abf0 --- /dev/null +++ b/examples/LArNEST/LegacyLArNESTBenchmarks.cpp @@ -0,0 +1,113 @@ +/** + * @file legacyLArNEST.cpp + * @author NEST Collaboration + * @author Nicholas Carrara [nmcarrara@ucdavis.edu] + * @brief + * @version + * @date 2022-04-14 + */ +#include +#include +#include +#include +#include + +#include "LArDetector.hh" +#include "LArNEST.hh" + +int main(int argc, char* argv[]) +{ + // this sample program takes several arguments + // including + // 1) num_events - the number of events to generate + // 2) pdgcode - the pdg code of the particle creating the deposition + // 3) density - the density of the LAr + // 4) track_length - the size of the track length for the deposition + // (optional) 5) seed - a seed for the random number generator + if (argc < 5) + { + std::cerr << "Error! This program expects six inputs!" << std::endl; + exit(0); + } + size_t num_events = atoi(argv[1]); + + // set up energy steps + int num_energy_steps = 50000; + std::vector energy_vals; + double start_val = .1; + double end_val = 1000; + double step_size = (end_val - start_val)/num_energy_steps; + for (size_t i = 0; i < num_energy_steps; i++) + { + energy_vals.emplace_back(start_val + step_size * i); + } + + int pdgcode = atoi(argv[2]); + double density = atof(argv[3]); + std::vector electric_field = { + 1, 50, 100, 200, 500, 1000, 1500, 2000 + }; + + double track_length = atof(argv[4]); + uint64_t seed = 0; + if (argc > 5) { + seed = atoi(argv[5]); + } + + LArDetector* detector = new LArDetector(); + // Construct NEST class using detector object + NEST::LArNEST larnest(detector); + // Set random seed of RandomGen class + RandomGen::rndm()->SetSeed(seed); + + // Construct NEST objects for storing calculation results + NEST::LArYieldResult result; + std::ofstream output_file; + output_file.open("legacy_larnest_benchmarks_" + std::to_string(pdgcode) + ".csv"); + output_file << "energy,efield,TotalYield,QuantaYield,LightYield,Nph,Ne,Nex,Nion\n"; + // iterate over number of events + for (size_t v = 0; v < electric_field.size(); v++) + { + for (size_t i = 0; i < num_energy_steps; i++) + { + std::vector TotalYield(num_events); + std::vector QuantaYield(num_events); + std::vector LightYield(num_events); + std::vector Nph(num_events); + std::vector Ne(num_events); + std::vector Nex(num_events); + std::vector Nion(num_events); + // collect statistics for each event + for (size_t j = 0; j < num_events; j++) + { + result = larnest.LegacyCalculation( + pdgcode, + energy_vals[i], + electric_field[v], + density, + track_length + ); + TotalYield[j] = result.TotalYield; + QuantaYield[j] = result.QuantaYield; + LightYield[j] = result.LightYield; + Nph[j] = result.Nph; + Ne[j] = result.Ne; + Nex[j] = result.Nex; + Nion[j] = result.Nion; + } + double TotalYield_mean = std::accumulate(TotalYield.begin(), TotalYield.end(), 0.0) / double(num_events); + double QuantaYield_mean = std::accumulate(QuantaYield.begin(), QuantaYield.end(), 0.0) / double(num_events); + double LightYield_mean = std::accumulate(LightYield.begin(), LightYield.end(), 0.0) / double(num_events); + double Nph_mean = std::accumulate(Nph.begin(), Nph.end(), 0.0) / double(num_events); + double Ne_mean = std::accumulate(Ne.begin(), Ne.end(), 0.0) / double(num_events); + double Nex_mean = std::accumulate(Nex.begin(), Nex.end(), 0.0) / double(num_events); + double Nion_mean = std::accumulate(Nion.begin(), Nion.end(), 0.0) / double(num_events); + output_file << energy_vals[i] << ","; + output_file << electric_field[v] << ","; + output_file << TotalYield_mean << "," << QuantaYield_mean << "," << LightYield_mean << ","; + output_file << Nph_mean << "," << Ne_mean << "," << Nex_mean << "," << Nion_mean << "\n"; + } + } + output_file.close(); + return 0; +} \ No newline at end of file diff --git a/include/Detectors/LArDetector.hh b/include/Detectors/LArDetector.hh new file mode 100644 index 00000000..1b4b4be4 --- /dev/null +++ b/include/Detectors/LArDetector.hh @@ -0,0 +1,24 @@ +/** + * @file LArDetector.hh + * @author NEST Collaboration + * @author Nicholas Carrara [nmcarrara@ucdavis.edu] + * @brief + * @version + * @date 2022-04-14 + */ +#pragma once + +#include "VDetector.hh" + +/** + * @brief An example LAr TPC detector. + * + */ +class LArDetector : public VDetector +{ +public: + LArDetector(); + ~LArDetector() override = default; + void Initialization() override; +private: +}; \ No newline at end of file diff --git a/include/NEST/LArNEST.hh b/include/NEST/LArNEST.hh new file mode 100644 index 00000000..748b476c --- /dev/null +++ b/include/NEST/LArNEST.hh @@ -0,0 +1,530 @@ +/** + * @file LArNEST.hh + * @author NEST Collaboration + * @author Nicholas Carrara [nmcarrara@ucdavis.edu] + * @author Justin Mueller [Justin.Mueller@colostate.edu] + * @author Ekaterina Kozlova [aspelene@gmail.com] + * @author Michael Mooney [mrmooney@colostate.edu] + * @brief + * @version + * @date 2022-04-13 + */ +#pragma once + +#include "NEST.hh" +#include "RandomGen.hh" + +namespace NEST +{ + static constexpr double LAr_Z{18}; + static constexpr double legacy_density_LAr{1.393}; + static constexpr double legacy_scint_yield{1.0 / (19.5 * 1.e-6)}; + static constexpr double legacy_resolution_scale{0.107}; // Doke 1976 + static constexpr double two_PI = 2. * M_PI; + static constexpr double sqrt2 = gcem::sqrt(2.); + static constexpr double sqrt2_PI = gcem::sqrt(2. * M_PI); + static constexpr double inv_sqrt2_PI = 1. / gcem::sqrt(2. * M_PI); + + enum class LArInteraction + { + NR = 0, + ER = 1, + Alpha = 2, + }; + + struct LArNRYieldsParameters + { + double alpha = {11.10}; + double beta = {0.087}; + double gamma = {0.1}; + double delta = {-0.0932}; + double epsilon ={2.998}; + double zeta = {0.3}; + double eta = {2.94}; + }; + struct LArERElectronYieldsAlphaParameters + { + double A = {32.988}; + double B = {-552.988}; + double C = {17.2346}; + double D = {-4.7}; + double E = {0.025115}; + double F = {0.265360653}; + double G = {0.242671}; + }; + struct LArERElectronYieldsBetaParameters + { + double A = {0.778482}; + double B = {25.9}; + double C = {1.105}; + double D = {0.4}; + double E = {4.55}; + double F = {-7.502}; + }; + struct LArERElectronYieldsGammaParameters + { + double A = {0.659509}; + double B = {1000}; + double C = {6.5}; + double D = {5.0}; + double E = {-0.5}; + double F = {1047.408}; + double G = {0.01851}; + }; + struct LArERElectronYieldsDokeBirksParameters + { + double A = {1052.264}; + double B = {14159350000 - 1652.264}; + double C = {-5.0}; + double D = {0.157933}; + double E = {1.83894}; + }; + struct LArERYieldsParameters + { + LArERElectronYieldsAlphaParameters alpha; + LArERElectronYieldsBetaParameters beta; + LArERElectronYieldsGammaParameters gamma; + LArERElectronYieldsDokeBirksParameters doke_birks; + double p1 = {1.0}; + double p2 = {10.304}; + double p3 = {13.0654}; + double p4 = {0.10535}; + double p5 = {0.7}; + double delta = {15.7489}; + double let = {-2.07763}; + }; + + struct LArAlphaElectronYieldsParameters + { + double A = {1.0/6200.0}; + double B = {64478398.7663}; + double C = {0.173553719}; + double D = {1.21}; + double E = {0.02852}; + double F = {0.01}; + double G = {4.71598}; + double H = {7.72848}; + double I = {-0.109802}; + double J = {3.0}; + }; + + struct LArAlphaPhotonYieldsParameters + { + double A = {1.5}; + double B = {-0.012}; + double C = {1.0/6500.0}; + double D = {278037.250283}; + double E = {0.173553719}; + double F = {1.21}; + double G = {2}; + double H = {0.653503}; + double I = {4.98483}; + double J = {10.0822}; + double K = {1.2076}; + double L = {-0.97977}; + double M = {3.0}; + }; + + struct LArAlphaYieldsParameters + { + LArAlphaElectronYieldsParameters Ye; + LArAlphaPhotonYieldsParameters Yph; + }; + + struct ThomasImelParameters + { + double A = {0.1}; + double B = {-0.0932}; + }; + + struct DriftParameters + { + std::vector A = { + 0.937729, 0.80302379, 0.7795972, 0.6911897, + 0.76551511, 0.502022794, 0.24207633 + }; + std::vector B = { + -0.0734108, -0.06694564, -0.0990952, -0.092997, + -0.0731659, -0.06644517, -0.03558428 + }; + std::vector C = { + 0.315338, 0.331798, 0.320876, 0.3295202, + 0.317972, 0.3290246, 0.33645519 + }; + + std::vector TempLow = { + 84., 86., 88., 92., 96., 110., 125. + }; + std::vector TempHigh = { + 86., 88., 92., 96., 110., 125., 140. + }; + }; + + struct LArYieldResult + { + double TotalYield; + double QuantaYield; + double LightYield; + double Nph; + double Ne; + double Nex; + double Nion; + double Lindhard; + double ElectricField; + }; + + struct LArYieldFluctuationResult + { + double TotalYieldFluctuation; + double QuantaYieldFluctuation; + double LightYieldFluctuation; + double NphFluctuation; + double NeFluctuation; + double NexFluctuation; + double NionFluctuation; + }; + + struct LArNESTResult + { + LArYieldResult yields; + LArYieldFluctuationResult fluctuations; + photonstream photon_times; + }; + + /** + * @brief + * + */ + class LArNEST : public NESTcalc + { + public: + explicit LArNEST(VDetector *detector); + + //-------------------------Parameters-------------------------// + /// set LAr parameters + void setDensity(double density) { fDensity = density; } + void setRIdealGas(double RIdealGas) { fRIdealGas = RIdealGas; } + void setRealGasA(double RealGasA) { fRealGasA = RealGasA; } + void setRealGasB(double RealGasB) { fRealGasB = RealGasB; } + void setWorkQuantaFunction(double workQuantaFunction) + { fWorkQuantaFunction = workQuantaFunction; } + void setWorkIonFunction(double workIonFunction) + { fWorkIonFunction = workIonFunction; } + void setWorkPhotonFunction(double workPhotonFunction) + { fWorkPhotonFunction = workPhotonFunction; } + void setFanoER(double FanoER) { fFanoER = FanoER; } + void setNexOverNion(double NexOverNion) { fNexOverNion = NexOverNion; } + + // TODO: do we need nuisance parameters? + void setNuisanceParameters(std::vector nuisanceParameters); + void setTemperature(std::vector temperature); + + /// setters for various parameters + void setNRYieldsParameters(LArNRYieldsParameters NRYieldsParameters); + void setERYieldsParameters(LArERYieldsParameters ERYieldsParameters); + void setERElectronYieldsAlphaParameters( + LArERElectronYieldsAlphaParameters ERElectronYieldsAlphaParameters + ); + void setERElectronYieldsBetaParameters( + LArERElectronYieldsBetaParameters ERElectronYieldsBetaParameters + ); + void setERElectronYieldsGammaParameters( + LArERElectronYieldsGammaParameters ERElectronYieldsGammaParameters + ); + void setERElectronYieldsDokeBirksParameters( + LArERElectronYieldsDokeBirksParameters ERElectronYieldsDokeBirksParameters + ); + void setThomasImelParameters(ThomasImelParameters thomasImelParameters); + void setDriftParameters(DriftParameters driftParameters); + + /// get LAr parameters + double getDensity() const { return fDensity; } + double getRIdealGas() const { return fRIdealGas; } + double getRealGasA() const { return fRealGasA; } + double getRealGasB() const { return fRealGasB; } + double getWorkQuantaFunction() const { return fWorkQuantaFunction; } + double getWorkIonFunction() const { return fWorkIonFunction; } + double getWorkPhotonFunction() const { return fWorkPhotonFunction; } + double getFanoER() const { return fFanoER; } + double getNexOverNion() const { return fNexOverNion; } + + LArNRYieldsParameters getNRYieldsParameters() { return fNR; } + LArERYieldsParameters getERYieldsParameters() { return fER; } + LArERElectronYieldsAlphaParameters getERElectronYieldsAlphaParameters() + { return fER.alpha; } + LArERElectronYieldsBetaParameters getERElectronYieldsBetaParameters() + { return fER.beta; } + LArERElectronYieldsGammaParameters getERElectronYieldsGammaParameters() + { return fER.gamma; } + LArERElectronYieldsDokeBirksParameters getERElectronYieldsDokeBirksParameters() + { return fER.doke_birks; } + ThomasImelParameters getThomasImelParameters() {return fThomasImelParameters; } + DriftParameters getDriftParameters() {return fDriftParameters; } + + //-------------------------All Yields-------------------------// + LArYieldResult GetRecombinationYields( + double TotalYields, double ElectronYields, double PhotonYields, + double energy, double efield + ); + LArYieldResult GetYields( + LArInteraction species, double energy, + double efield, double density + ); + /** + * @brief Calculate fluctions on the mean yields + * + */ + LArYieldFluctuationResult GetYieldFluctuations( + const YieldResult &yields, double density + ); + /** + * @brief + * + */ + LArNESTResult FullCalculation( + LArInteraction species, double energy, + double efield, double density, + bool do_times + ); + + //-------------------------NR Yields-------------------------// + /** + * @brief NR Total Yields function: + * Ty = alpha * E^beta + * + * @param energy + * @return double + */ + double GetNRTotalYields(double energy); + /** + * @brief NR Exciton Yields function: + * Qy = (gamma * E^delta)*(sqrt(E + eps)^-1) + * * (1 - (1 + (E/zeta)^eta))^-1) + * + * @param energy + * @param efield + * @return double + */ + double GetNRElectronYields(double energy, double efield); + /** + * @brief NR Photon Yields function: + * Ly = + * + * @param energy + * @param efield + * @return double + */ + double GetNRPhotonYields(double energy, double efield); + /** + * @brief NR Photon Yields function (conserved) + * Ly = Ty/E - Qy + * + * @param energy + * @param efield + * @return double + */ + double GetNRPhotonYieldsConserved(double energy, double efield); + /** + * @brief Calculate yields for nuclear recoils. The formulas + * for these are given by: + * Ty = alpha * E^beta + * Qy = (gamma * F^delta)^-1 * (sqrt(E + epsilon)) + * * (1 - (1 + (E/zeta)^eta)^-1) + * + * Ly = alpha * E^(beta - 1) - (gamma * F^delta)^-1 + * * (sqrt(E + epsilon))^-1 + * First, we fit the total yield model and take it as fixed. + * The model is a simple power law in the deposited energy. There + * is no field dependence. (J. Mueller and E. Kozlova) + * @param energy + * @param density + * @param efield + * @param NuisParam + * @return YieldResult + */ + LArYieldResult GetNRYields( + double energy, double efield, double density + ); + + //-------------------------ER Yields-------------------------// + /** + * @brief ER Total Yields function: + * Ty + * + * @param energy + * @return double + */ + double GetERTotalYields(double energy); + /** + * @brief ER Exciton Yields Alpha function: + * + * @param efield + * @param density + * @return double + */ + double GetERElectronYieldsAlpha(double efield, double density); + /** + * @brief ER Exciton Yields Beta function: + * + * @param efield + * @return double + */ + double GetERElectronYieldsBeta(double efield); + /** + * @brief ER Exciton Yields Gamma function: + * + * @param efield + * @return double + */ + double GetERElectronYieldsGamma(double efield); + /** + * @brief ER Exciton Yields Doke-Birks function: + * + * @param efield + * @return double + */ + double GetERElectronYieldsDokeBirks(double efield); + /** + * @brief ER Exciton Yields function: + * + * @param energy + * @param efield + * @param density + * @return double + */ + double GetERElectronYields(double energy, double efield, double density); + /** + * @brief Calculate yields for charged particles. The formulas + * for these are given by: + * Qy = alpha * beta + (gamma - alpha * beta)/(p1 + p2 * E^p3)^p4 + * + delta / (p5 + Doke * E^LET) + * + * Ly = Nq - Qy + * + * The p1-p5 parameters are stored in ERQuantaParameters. + * @param energy + * @param density + * @param efield + * @return LArYieldResult + */ + LArYieldResult GetERYields( + double energy, double efield, double density + ); + //-------------------------Alpha Yields-------------------------// + /** + * @brief Get the Alpha Total Yields object + * + * @param energy + * @return double + */ + double GetAlphaTotalYields(double energy); + /** + * @brief Get the Alpha Electron Yields object + * + * @param efield + * @return double + */ + double GetAlphaElectronYields(double efield); + /** + * @brief Get the Alpha Photon Yields object + * + * @param efield + * @return double + */ + double GetAlphaPhotonYields(double efield); + /** + * @brief Calculate yields for alpha particles + * + * @param energy + * @param efield + * @param density + * @return LArYieldResult + */ + LArYieldResult GetAlphaYields( + double energy, double efield, double density + ); + + //-------------------------Photon Times-------------------------// + double GetPhotonTime( + LArInteraction species, bool exciton, + double energy + ); + double GetPhotonEnergy(bool state); + + //-------------------------Drift Velocity-------------------------// + double GetDriftVelocity_Liquid( + double Kelvin, double eField + ); + double GetDriftVelocity_MagBoltz( + double density, double efieldinput, + double molarMass + ); + + //-------------------------Utilities-------------------------// + double GetLinearEnergyTransfer( + double energy, bool CSDA=false + ); + double GetDensity( + double Kelvin, double bara, bool &inGas, + uint64_t evtNum, double molarMass + ); + std::vector CalculateG2(bool verbosity=false); + + //-------------------------Legacy LArNEST-------------------------// + /** + * @brief + * + * @param energy + * @param efield + * @param yieldFactor + * @param excitationRatio + * @param epsilon + * @param recombProb + * @return LArYieldResult + */ + LArYieldResult LegacyGetYields( + double energy, double efield, + double yieldFactor, double excitationRatio, + double epsilon, double recombProb + ); + /** + * @brief Below are legacy LAr calculation + * functions which are almost copied verbadim + * from LArSoft, see - + * + * The default value for the dx=track_length is + * 1/3 mm, which is the standard from Geant4. + * @return LArYieldResult + */ + LArYieldResult LegacyCalculation( + int pdgcode, double energy, + double efield, double density, + double track_length=0.0003 + ); + double LegacyGetRecombinationProbability( + double energy, double efield, double density, + int pdgcode, double track_length + ); + double LegacyGetLinearEnergyTransfer(double E); + + private: + double fDensity = {1.393}; + double fRIdealGas = {8.31446261815324}; + double fRealGasA = {0.1355}; // m^6*Pa/mol^2 or m^4*N/mol^2. + double fRealGasB = {3.201e-5}; // m^3/mol. + + double fWorkQuantaFunction = {19.5}; + double fWorkIonFunction = {23.6}; + double fWorkPhotonFunction = {14.544}; + + double fNexOverNion = {0.21}; + double fFanoER = {0.1115}; + + LArNRYieldsParameters fNR; + LArERYieldsParameters fER; + LArAlphaYieldsParameters fAlpha; + + ThomasImelParameters fThomasImelParameters; + DriftParameters fDriftParameters; + }; +} \ No newline at end of file diff --git a/include/NEST/TestLArSpectra.hh b/include/NEST/TestLArSpectra.hh new file mode 100644 index 00000000..6e791b0e --- /dev/null +++ b/include/NEST/TestLArSpectra.hh @@ -0,0 +1,54 @@ +/** + * @file TestLArSpectra.hh + * @author NEST Collaboration + * @author Nicholas Carrara [nmcarrara@ucdavis.edu] + * @brief + * @version + * @date 2022-04-14 + */ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "ValidityTests.hh" +#include "RandomGen.hh" + +class TestLArSpectra { + public: + TestLArSpectra() = default; + // struct WIMP_spectrum_prep { + // double base[100] = {1.}; + // double exponent[100] = {0.}; + // double integral = 0.; + // double xMax = 0.; + // double divisor = 1.; + // }; + // WIMP_spectrum_prep wimp_spectrum_prep; + + // static double CH3T_spectrum(double emin, double emax); + // static double C14_spectrum(double emin, double emax); + // static double B8_spectrum(double emin, double emax); + // static double AmBe_spectrum(double emin, double emax); + // static double Cf_spectrum(double emin, double emax); + // static double DD_spectrum(double xMin, double xMax, double expFall, double peakFrac, double peakMu, double peakSig); + // static double ppSolar_spectrum(double emin, double emax); + // static double atmNu_spectrum(double emin, double emax); + // static double WIMP_dRate(double ER, double mWimp, double day); + // static WIMP_spectrum_prep WIMP_prep_spectrum(double mass, double eStep, + // double day); + // static double WIMP_spectrum(WIMP_spectrum_prep wprep, double mass, + // double day); + // static const vector Gamma_spectrum(double xMin, double xMax, + // string source); + + // static double + // ZeplinBackground(); // an example of how to do a better (non-flat) ER + // // BG spectrum for a WS, from Henrique Araujo +}; \ No newline at end of file diff --git a/src/LArDetector.cpp b/src/LArDetector.cpp new file mode 100644 index 00000000..585cdfaf --- /dev/null +++ b/src/LArDetector.cpp @@ -0,0 +1,87 @@ +/** + * @file LArDetector.cpp + * @author NEST Collaboration + * @author Nicholas Carrara [nmcarrara@ucdavis.edu] + * @brief + * @version + * @date 2022-04-14 + */ +#include "LArDetector.hh" + +LArDetector::LArDetector() +{ + Initialization(); +} + +void LArDetector::Initialization() +{ + // Primary Scintillation (S1) parameters + // phd per S1 phot at dtCntr (not phe). Divide out 2-PE effect + double g1 = 0.0760; + // single phe resolution (Gaussian assumed) + double sPEres = 0.58; + // POD threshold in phe, usually used IN PLACE of sPEeff + double sPEthr = 0.35; + // actual efficiency, can be used in lieu of POD threshold + double sPEeff = 1.00; + double noiseBaseline[4] = { + 0.0, // baseline noise mean and width in PE (Gaussian) + 0.0, // baseline noise mean and width in PE (Gaussian) + 0.0, // EXO noise mean + 0.0 // EXO noise width + }; + // chance 1 photon makes 2 phe instead of 1 in Hamamatsu PMT + double P_dphe = 0.2; + + double coinWind = 100; // S1 coincidence window in ns + int coinLevel = 2; // how many PMTs have to fire for an S1 to count + int numPMTs = 89; // For coincidence calculation + + bool OldW13eV = true; // for matching EXO-200's W measurement + // "Linear noise" terms as defined in Dahl thesis and by D. McK + // S1->S1 Gaussian-smeared w/ noiseL[0]*S1. Ditto S2 + double noiseLinear[2] = { + 3e-2, 3e-2 + }; + //(n)EXO quadratic noise term + double noiseQuadratic[2] = {0.0, 0.0}; + + // Ionization and Secondary Scintillation (S2) parameters + double g1_gas = 0.06; // phd per S2 photon in gas, used to get SE size + double s2Fano = 3.61; // Fano-like fudge factor for SE width + // the S2 threshold in phe or PE, *not* phd. Affects NR most + double s2_thr = 300.; + double E_gas = 12.; // field in kV/cm between liquid/gas border and anode + double eLife_us = 2200.; // the drift electron mean lifetime in micro-seconds + + // Thermodynamic Properties + bool inGas = false; + double T_Kelvin = 177.; // for liquid drift speed calculation + double p_bar = 2.14; // gas pressure in units of bars, it controls S2 size + // if you are getting warnings about being in gas, lower T and/or raise p + + // Data Analysis Parameters and Geometry + double dtCntr = 40.; // center of detector for S1 corrections, in usec. + double dt_min = 20.; // minimum. Top of detector fiducial volume + double dt_max = 60.; // maximum. Bottom of detector fiducial volume + + double radius = 50.; // millimeters + double radmax = 50.; + + double TopDrift = 150.; // mm not cm or us (but, this *is* where dt=0) + // a z-axis value of 0 means the bottom of the detector (cathode OR bottom + // PMTs) + // In 2-phase, TopDrift=liquid/gas border. In gas detector it's GATE, not + // anode! + double anode = 152.5; // the level of the anode grid-wire plane in mm + // In a gas TPC, this is not TopDrift (top of drift region), but a few mm + // above it + double gate = 147.5; // mm. This is where the E-field changes (higher) + // in gas detectors, the gate is still the gate, but it's where S2 starts + double cathode = 1.00; // mm. Defines point below which events are gamma-X + + // 2-D (X & Y) Position Reconstruction + double PosResExp = 0.015; // exp increase in pos recon res at hi r, 1/mm + double PosResBase = 70.8364; // baseline unc in mm, see NEST.cpp for usage + double molarMass = 131.293; // molar mass, g/mol +} \ No newline at end of file diff --git a/src/LArNEST.cpp b/src/LArNEST.cpp new file mode 100644 index 00000000..ef9c8170 --- /dev/null +++ b/src/LArNEST.cpp @@ -0,0 +1,612 @@ +/** + * @file LArNEST.cpp + * @author NEST Collaboration + * @author Nicholas Carrara [nmcarrara@ucdavis.edu] + * @author Justin Mueller [Justin.Mueller@colostate.edu] + * @author Ekaterina Kozlova [aspelene@gmail.com] + * @author Michael Mooney [mrmooney@colostate.edu] + * @brief + * @version + * @date 2022-04-13 + */ +#include "LArNEST.hh" + +namespace NEST +{ + LArNEST::LArNEST(VDetector *detector) + : NESTcalc(detector) + {} + //-------------------------All Yields-------------------------// + LArYieldResult LArNEST::GetRecombinationYields( + double TotalYields, double ElectronYields, double PhotonYields, + double energy, double efield + ) + { + double Ne = ElectronYields * energy; + double Nph = PhotonYields * energy; + double Nq = Ne + Nph; + + /// recombination + double ThomasImel = fThomasImelParameters.A * pow(efield, fThomasImelParameters.B); + double Nion = (4. / ThomasImel) * (exp(Ne * ThomasImel / 4.) - 1.); + if (Nion < 0.0) { + Nion = 0.0; + } + if (Nion > Nq) { + Nion = Nq; + } + double Nex = Nq - Nion; + + double WQ_eV = fWorkQuantaFunction; + double Lindhard = (TotalYields / energy) * WQ_eV * 1e-3; + LArYieldResult result{ + TotalYields, ElectronYields, PhotonYields, + Nph, Ne, Nex, Nion, Lindhard, efield + }; + return result; + } + LArYieldResult LArNEST::GetYields( + LArInteraction species, double energy, + double efield, double density + ) + { + if (species == LArInteraction::NR) { + return GetNRYields(energy, efield, density); + } + else if (species == LArInteraction::ER) { + return GetERYields(energy, efield, density); + } + else { + return GetAlphaYields(energy, efield, density); + } + } + + LArYieldFluctuationResult LArNEST::GetYieldFluctuations( + const YieldResult &yields, double density + ) + { + // TODO: Understand and break up this function. + LArYieldFluctuationResult result{}; + return result; // quanta returned with recomb fluctuations + } + + LArNESTResult LArNEST::FullCalculation( + LArInteraction species, double energy, + double efield, double density, + bool do_times + ) + { + if (density < 1.) fdetector->set_inGas(true); + LArNESTResult result; + result.yields = GetYields(species, energy, efield, density); + // result.quanta = GetQuanta(result.yields, density); + // if (do_times) + // result.photon_times = GetPhotonTimes( + // species, result.quanta.photons, + // result.quanta.excitons, efield, energy + // ); + // else { + // result.photon_times = photonstream(result.quanta.photons, 0.0); + // } + return result; + } + + //-------------------------NR Yields-------------------------// + double LArNEST::GetNRTotalYields(double energy) + { + return fNR.alpha * pow(energy, fNR.beta); + } + double LArNEST::GetNRElectronYields(double energy, double efield) + { + return (1.0 / (fNR.gamma * pow(efield, fNR.delta))) * + (1.0 / sqrt(energy + fNR.epsilon)) * + (1.0 - 1.0 / (1.0 + pow(energy/fNR.zeta, fNR.eta))); + } + double LArNEST::GetNRPhotonYields(double energy, double efield) + { + return fNR.alpha * pow(energy, fNR.beta) - + (1.0 / (fNR.gamma * pow(efield, fNR.delta))) * + (1.0 / sqrt(energy + fNR.epsilon)); + } + double LArNEST::GetNRPhotonYieldsConserved(double energy, double efield) + { + return fNR.alpha * pow(energy, fNR.beta) - + (1.0 / (fNR.gamma * pow(efield, fNR.delta))) * + (1.0 / sqrt(energy + fNR.epsilon)) * + (1.0 - 1.0 / (1.0 + pow(energy/fNR.zeta, fNR.eta))); + } + LArYieldResult LArNEST::GetNRYields( + double energy, double efield, double density + ) + { + double NRTotalYields = GetNRTotalYields(energy); + double NRElectronYields = GetNRElectronYields(energy, efield); + if (NRElectronYields < 0.0) { + NRElectronYields = 0.0; + } + double NRPhotonYields = GetNRPhotonYields(energy, efield); + if (NRPhotonYields < 0.0) { + NRPhotonYields = 0.0; + } + return GetRecombinationYields( + NRTotalYields, NRElectronYields, NRPhotonYields, + energy, efield + ); + } + + //-------------------------ER Yields-------------------------// + double LArNEST::GetERTotalYields(double energy) + { + return (energy * 1.0e3 / fWorkQuantaFunction); + } + double LArNEST::GetERElectronYieldsAlpha(double efield, double density) + { + return fER.alpha.A + + fER.alpha.B / + (fER.alpha.C + + pow(efield / (fER.alpha.D + fER.alpha.E * + exp(density / fER.alpha.F)), fER.alpha.G)); + } + double LArNEST::GetERElectronYieldsBeta(double efield) + { + return fER.beta.A + + fER.beta.B * + pow(fER.beta.C + + pow(efield / fER.beta.D, fER.beta.E), fER.beta.F); + } + double LArNEST::GetERElectronYieldsGamma(double efield) + { + return fER.gamma.A * + (fER.gamma.B / fWorkQuantaFunction + + fER.gamma.C * + (fER.gamma.D + fER.gamma.E / + pow(efield / fER.gamma.F, fER.gamma.G))); + } + double LArNEST::GetERElectronYieldsDokeBirks(double efield) + { + return fER.doke_birks.A + + fER.doke_birks.B / + (fER.doke_birks.C + + pow(efield / fER.doke_birks.D, fER.doke_birks.E)); + } + double LArNEST::GetERElectronYields(double energy, double efield, double density) + { + double alpha = GetERElectronYieldsAlpha(efield, density); + double beta = GetERElectronYieldsBeta(efield); + double gamma = GetERElectronYieldsGamma(efield); + double doke_birks = GetERElectronYieldsDokeBirks(efield); + return alpha * beta + + (gamma - alpha * beta) / + pow(fER.p1 + fER.p2 * pow(energy + 0.5, fER.p3), fER.p4) + + fER.delta / (fER.p5 + doke_birks * pow(energy, fER.let)); + } + LArYieldResult LArNEST::GetERYields( + double energy, double efield, double density + ) + { + double ERTotalYields = GetERTotalYields(energy); + double ERElectronYields = GetERElectronYields(energy, efield, density); + if (ERElectronYields < 0.0) { + ERElectronYields = 0.0; + } + double ERPhotonYields = ERTotalYields / energy - ERElectronYields; + if (ERPhotonYields < 0.0) { + ERPhotonYields = 0.0; + } + return GetRecombinationYields( + ERTotalYields, ERElectronYields, ERPhotonYields, + energy, efield + ); + } + + //-------------------------Alpha Yields-------------------------// + double LArNEST::GetAlphaTotalYields(double energy) + { + } + double LArNEST::GetAlphaElectronYields(double efield) + { + double fieldTerm = fAlpha.Ye.F * pow( + (fAlpha.Ye.G + pow(efield, fAlpha.Ye.H)),fAlpha.Ye.I + ); + return fAlpha.Ye.A * ( + fAlpha.Ye.B - (fAlpha.Ye.B * fAlpha.Ye.C + (fAlpha.Ye.B / fAlpha.Ye.D) * ( + 1 - (fAlpha.Ye.E * log(1 + (fAlpha.Ye.B / fAlpha.Ye.D) * (fieldTerm / fAlpha.Ye.J)) / ( + (fAlpha.Ye.B / fAlpha.Ye.D) * fieldTerm) + ))) + ); + } + double LArNEST::GetAlphaPhotonYields(double efield) + { + double quench = 1.0 / (fAlpha.Yph.A * pow(efield, fAlpha.Yph.B)); + double fieldTerm = fAlpha.Yph.H * pow( + (fAlpha.Yph.I + pow((efield/fAlpha.Yph.J), fAlpha.Yph.K)),fAlpha.Yph.L + ); + return quench * fAlpha.Yph.C * ( + fAlpha.Yph.D * fAlpha.Yph.E + (fAlpha.Yph.D / fAlpha.Yph.F) * ( + 1 - (fAlpha.Yph.G * log(1 + (fAlpha.Yph.D / fAlpha.Yph.F) * (fieldTerm / fAlpha.Yph.M)) / ( + (fAlpha.Yph.D / fAlpha.Yph.F) * fieldTerm) + )) + ); + } + LArYieldResult LArNEST::GetAlphaYields( + double energy, double efield, double density + ) + { + double AlphaElectronYields = GetAlphaElectronYields(efield); + if (AlphaElectronYields < 0.0) { + AlphaElectronYields = 0.0; + } + double AlphaPhotonYields = GetAlphaPhotonYields(efield); + if (AlphaPhotonYields < 0.0) { + AlphaPhotonYields = 0.0; + } + double AlphaTotalYields = (AlphaElectronYields + AlphaPhotonYields) * energy; + return GetRecombinationYields( + AlphaTotalYields, AlphaElectronYields, AlphaPhotonYields, + energy, efield + ); + } + + //-------------------------Photon Times-------------------------// + double LArNEST::GetPhotonTime( + LArInteraction species, bool exciton, + double energy + ) + { + // arXiv:1802.06162. NR may need tauR ~0.5-1ns instead of 0 + double SingTripRatio = 0.; + // if ( efield > 60. ) efield = 60. // makes Xed work. 200 for LUX Run04 + // instead. A mystery! Why no field dep? + + double LET = GetLinearEnergyTransfer(energy); + // copied from 2013 NEST version for LAr on LBNE + double tau1 = RandomGen::rndm()->rand_gauss(6.5, 0.8); // error from weighted average + double tau3 = RandomGen::rndm()->rand_gauss(1300, 50); // ibid. + double tauR = RandomGen::rndm()->rand_gauss(0.8, 0.2); // Kubota 1979 + if (energy < fWorkQuantaFunction * 0.001) + { + // from old G4S2Light + tau1 = 6.; + tau3 = 1600.; + SingTripRatio = 0.1; + if (!exciton) { + tauR = 28e3; + } + else { + tauR = 0.; // 28 microseconds comes from Henrique: + // https://doi.org/10.1016/j.astropartphys.2018.04.006 + } + } + else + { + if (species == LArInteraction::NR) { + SingTripRatio = 0.22218 * pow(energy, 0.48211); + } + else if (species == LArInteraction::Alpha) + { // really only alphas here + SingTripRatio = (-0.065492 + 1.9996 * exp(-energy / 1e3)) / + (1. + 0.082154 / pow(energy / 1e3, 2.)) + + 2.1811; // uses energy in MeV not keV + } + else + { + if (LET < 3.) + { + if (!exciton) { + SingTripRatio = RandomGen::rndm()->rand_gauss(0.5, 0.2); + } + else if (exciton) { + SingTripRatio = RandomGen::rndm()->rand_gauss(0.36, 0.06); + } + } + else + { + SingTripRatio = 0.2701 + + 0.003379 * LET - + 4.7338e-5 * pow(LET, 2.) + + 8.1449e-6 * pow(LET, 3.); + } + } // lastly is ER for LAr + } + // in case varied with Gaussian earlier + // the recombination time is non-exponential, but approximates + // to exp at long timescales (see Kubota '79) + double time_ns = tauR * (1.0 / RandomGen::rndm()->rand_uniform() - 1.); + + if (RandomGen::rndm()->rand_uniform() < SingTripRatio / (1. + SingTripRatio)) { + time_ns -= tau1 * log(RandomGen::rndm()->rand_uniform()); + } + else { + time_ns -= tau3 * log(RandomGen::rndm()->rand_uniform()); + } + return time_ns; + } + double LArNEST::GetPhotonEnergy(bool state) + { + // liquid Argon + // TODO: What is this function about? + if (state) { + return RandomGen::rndm()->rand_gauss(9.7, 0.2); + } + else { + return RandomGen::rndm()->rand_gauss(9.69, 0.22); + } + } + + //-------------------------Drift Velocity-------------------------// + double LArNEST::GetDriftVelocity_Liquid( + double Kelvin, double efield + ) + { + // returns drift speed in mm/usec. based on Fig. 14 arXiv:1712.08607 + if (Kelvin < 84. || Kelvin > 140.) + { + cerr << "\nWARNING: TEMPERATURE OUT OF RANGE (84-140 K) for vD\n"; + cerr << "Using value at closest temp for a drift speed estimate\n"; + Kelvin = (double)NESTcalc::clamp(int(Kelvin), 84, 140); + } + for (size_t i = 0; i < fDriftParameters.A.size() - 1; i++) + { + if (Kelvin >= fDriftParameters.TempLow[i] and Kelvin < fDriftParameters.TempHigh[i+1]) + { + double speed = exp( + fDriftParameters.A[i] + fDriftParameters.B[i] / (efield / 1000) + + fDriftParameters.C[i] * log(efield / 1000) + ); + if (speed > 0.0) { + return speed; + } + else { + return 0.0; + } + } + } + return 0.0; + } + + //-------------------------Utilities-------------------------// + double LArNEST::GetDensity( + double Kelvin, double bara, bool &inGas, + uint64_t evtNum, double molarMass + ) + { + // currently only for fixed pressure + // (the saturated vapor pressure); will add + // pressure dependence later + double VaporP_bar = exp( + 45.973940 - (1464.718291 / Kelvin) - (6.539938 * log(Kelvin)) + ); + if (bara < VaporP_bar || inGas) + { + double p_Pa = bara * 1e5; + double density = + 1. / + (pow(fRIdealGas * Kelvin, 3.) / + (p_Pa * pow(fRIdealGas * Kelvin, 2.) + fRealGasA * p_Pa * p_Pa) + + fRealGasB) * molarMass * 1e-6; // Van der Waals equation, mol/m^3 + if (bara < VaporP_bar && evtNum == 0) + // TODO: get rid of these output statements + std::cerr << "\nWARNING: ARGON GAS PHASE. IS THAT WHAT YOU WANTED?\n"; + inGas = true; + return density; + } + else + { + inGas = false; + return 1.4; + } + } + double LArNEST::GetLinearEnergyTransfer(double energy, bool CSDA) + { + double LET; + if (!CSDA) + { //total stopping power directly from ESTAR (radiative + collision) + LET = 1.8106 - + 0.45086 * log10(energy) - + 0.33151 * pow(log10(energy),2.) + + 0.25916 * pow(log10(energy),3.) - + 0.2051 * pow(log10(energy),4.) + + 0.15279 * pow(log10(energy),5.) - + 0.084659 * pow(log10(energy),6.) + + 0.030441 * pow(log10(energy),7.) - + 0.0058953 * pow(log10(energy),8.) + + 0.00045633 * pow(log10(energy),9.); + LET = pow(10.,LET); + if ( std::isnan(LET) || LET <= 0. ) { + LET = 1e2; + } + } + else + { + // the "continuous slowing down approximation" (CSDA) + // replace with Justin and Prof. Mooney's work + if (energy >= 1.) + { + LET = 116.70 - + 162.97 * log10(energy) + + 99.361 * pow(log10(energy), 2) - + 33.405 * pow(log10(energy), 3) + + 6.5069 * pow(log10(energy), 4) - + 0.69334 * pow(log10(energy), 5) + + 0.031563 * pow(log10(energy), 6); + } + else if (energy > 0. && energy < 1.) { + LET = 100.; + } + else { + LET = 0.; + } + } + return LET; + } + std::vector LArNEST::CalculateG2(bool verbosity) + { + } + + //------------------------------------Legacy LArNEST------------------------------------// + LArYieldResult LArNEST::LegacyGetYields( + double energy, double efield, + double yieldFactor, double excitationRatio, + double epsilon, double recombProb + ) + { + // determine ultimate number of quanta from current E-deposition (ph+e-) + // total mean number of exc/ions the total number of either quanta produced + // is equal to product of the work function, the energy deposited, + // and yield reduction, for NR + double MeanNq = legacy_scint_yield * energy; + double sigma = sqrt(legacy_resolution_scale * MeanNq); //Fano + double leftvar = RandomGen::rndm()->rand_gauss(yieldFactor, 0.25 * yieldFactor); + if (leftvar < 0) { + leftvar = 0; + } + if (leftvar > 1.0) { + leftvar = 1.0; + } + int Nq = int(floor(RandomGen::rndm()->rand_gauss(MeanNq, sigma) + 0.5)); + if (yieldFactor < 1) { + Nq = RandomGen::rndm()->binom_draw(Nq, leftvar); + } + // if Edep below work function, can't make any quanta, and if Nq + // less than zero because Gaussian fluctuated low, update to zero + if (energy < 1 / legacy_scint_yield || Nq < 0) { + Nq = 0; + } + + // next section binomially assigns quanta to excitons and ions + double Nex = RandomGen::rndm()->binom_draw( + Nq, excitationRatio / (1 + excitationRatio) + ); + double Nion = Nq - Nex; + + //use binomial distribution to assign photons, electrons, where photons + //are excitons plus recombined ionization electrons, while final + //collected electrons are the "escape" (non-recombined) electrons + double Nph = Nex + RandomGen::rndm()->binom_draw(Nion, recombProb); + double Ne = Nq - Nph; + + // create the quanta results + LArYieldResult result { + (Ne + Nph)/energy, Ne/energy, Nph/energy, + Nph, Ne, Nex, Nion, 0.0, efield + }; + return result; + } + LArYieldResult LArNEST::LegacyCalculation( + int pdgcode, double energy, + double efield, double density, + double track_length + ) + { + // determine various parameter values, such as the + // yieldfactor, excitationratio, dokebirks parameters, + // epislon and the recombination probability + // default quenching factor, for electronic recoils + double yieldFactor = 1.0; + // ratio for light particle in LAr, such as e-, mu-, Aprile et. al book + double excitationRatio = 0.21; + + // nuclear recoil quenching "L" factor: total yield is + // reduced for nuclear recoil as per Lindhard theory + double epsilon = 11.5 * (energy * pow(LAr_Z, (-7. / 3.))); + if (abs(pdgcode) == 2112) //nuclear recoil + { + yieldFactor = 0.23 * (1 + exp(-5 * epsilon)); //liquid argon L_eff + excitationRatio = 0.69337 + 0.3065 * exp(-0.008806 * pow(efield, 0.76313)); + } + // get the recombination probability + double recombProb = LegacyGetRecombinationProbability( + energy, efield, density, pdgcode, track_length + ); + return LegacyGetYields( + energy, efield, + yieldFactor, excitationRatio, + epsilon, recombProb + ); + } + double LArNEST::LegacyGetRecombinationProbability( + double energy, double efield, double density, + int pdgcode, double track_length + ) + { + // this section calculates recombination following the modified + // Birks'Law of Doke, deposition by deposition, may be overridden + // later in code if a low enough energy necessitates switching to the + // Thomas-Imel box model for recombination instead (determined by site) + double dE = energy; + double dx = 0.0; + double LET = 0.0; + double recombProb; + if (abs(pdgcode) != 11 && abs(pdgcode) != 13) + { + // e-: 11, e+: -11, mu-: 13, mu+: -13 + // in other words, if it's a gamma,ion,proton,alpha,pion,et al. do not + // use the step length provided by Geant4 because it's not relevant, + // instead calculate an estimated LET and range of the electrons that + // would have been produced if Geant4 could track them + LET = LegacyGetLinearEnergyTransfer(dE); + if (LET) { + //find the range based on the LET + dx = (dE / 1e3) / (density * LET); + } + if (abs(pdgcode) == 2112) // nuclear recoils + { + dx = 0; + } + } + else //normal case of an e-/+ energy deposition recorded by Geant + { + dx = track_length / 10.; + if (dx) + { + LET = ((dE / 1e3) / dx) * (1 / density); //lin. energy xfer (prop. to dE/dx) + } + if (LET > 0 && dE > 0 && dx > 0) + { + double ratio = LegacyGetLinearEnergyTransfer(dE) / LET; + if (ratio < 0.7 && pdgcode == 11) + { + dx /= ratio; + LET *= ratio; + } + } + } + // set up DokeBirks coefficients + double DokeBirksA = 0.07 * pow((efield / 1.0e3), -0.85); + double DokeBirksC = 0.00; + if (efield == 0.0) { + DokeBirksA = 0.0003; + DokeBirksC = 0.75; + } + //B=A/(1-C) (see paper) + double DokeBirksB = DokeBirksA / (1 - DokeBirksC); + recombProb = ((DokeBirksA * LET) / (1 + DokeBirksB * LET) + + DokeBirksC) * (density / legacy_density_LAr); + + //check against unphysicality resulting from rounding errors + return clamp(recombProb, 0., 1.); + } + double LArNEST::LegacyGetLinearEnergyTransfer(double E) + { + double LET; + if (E >= 1) + { + LET = 116.70 - + 162.97 * log10(E) + + 99.361 * pow(log10(E), 2) - + 33.405 * pow(log10(E), 3) + + 6.5069 * pow(log10(E), 4) - + 0.69334 * pow(log10(E), 5) + + 0.031563 * pow(log10(E), 6); + } + else if (E > 0 && E < 1) + { + LET = 100; + } + else + { + LET = 0; + } + return LET; + } +} \ No newline at end of file diff --git a/src/NEST.cpp b/src/NEST.cpp index dbd1793b..fcb57ef9 100644 --- a/src/NEST.cpp +++ b/src/NEST.cpp @@ -1164,8 +1164,10 @@ YieldResult NESTcalc::GetYields( YieldResult NESTcalc::YieldResultValidity(YieldResult &res, const double energy, const double Wq_eV) { - assert(res.ElectronYield != -999 && res.PhotonYield != -999 && - res.ExcitonRatio != -999); + assert(res.ElectronYield != -999 && + res.PhotonYield != -999 && + res.ExcitonRatio != -999 + ); if (res.PhotonYield > energy / W_SCINT) res.PhotonYield = energy / W_SCINT; // yields can never exceed 1 / [ W ~ few eV ] diff --git a/src/TestLArSpectra.cpp b/src/TestLArSpectra.cpp new file mode 100644 index 00000000..c6399c24 --- /dev/null +++ b/src/TestLArSpectra.cpp @@ -0,0 +1,10 @@ +/** + * @file TestLArSpectra.cpp + * @author NEST Collaboration + * @author Nicholas Carrara [nmcarrara@ucdavis.edu] + * @brief + * @version + * @date 2022-04-14 + */ +#include "TestLArSpectra.hh" + diff --git a/utils/larnest/lar_dataset.py b/utils/larnest/lar_dataset.py new file mode 100644 index 00000000..b5cf19cd --- /dev/null +++ b/utils/larnest/lar_dataset.py @@ -0,0 +1,255 @@ +""" +Functions for constructing training datasets from +LAr data files. +""" +import numpy as np +import pandas as pd +from matplotlib import pyplot as plt +import csv +import os + +from parameters import * +from utils import generate_plot_grid + +class LArDataset: + """ + """ + def __init__(self, + dataset_file: str='data/lar_data.npz', + benchmarks_file: str='', + plot_config: dict=default_plot_config, + ): + self.dataset_file = dataset_file + self.benchmarks_file = benchmarks_file + self.plot_config = plot_config + # set up datasets + self.data = dict(np.load(self.dataset_file, allow_pickle=True)) + self.datasets = {key: self.data[key].item() for key in self.data} + + self.benchmarks = pd.read_csv( + self.benchmarks_file, + names=[ + 'type','energy', 'efield' , + 'TotalYields', 'QuantaYields', 'LightYields', + 'Nph' , 'Ne', 'Nex', 'Nion' + ], + header=0, + #dtype=float + ) + self.benchmarks.replace([np.inf, -np.inf, np.nan], 0.0) + self.ylabels = { + 'nr_total': r"$Y_q$ [quanta/keV]", + 'nr_charge': r"$Y_{e^-}$ [electron/keV]", + 'nr_light': r"$Y_{\gamma}$ [photon/keV]", + 'nr_nph': r"$\langle N_{ph}\rangle$ [photon]", + 'nr_ne': r"$\langle N_{e^{-}}\rangle$ [electron]", + 'nr_nex': r"$\langle N_{ex}\rangle$ [exciton]", + 'nr_nion': r"$\langle N_{ion}\rangle$ [ion]", + 'er_total': r"$Y_q$ [quanta/keV]", + 'er_charge': r"$Y_{e^-}$ [electron/keV]", + 'er_light': r"$Y_{\gamma}$ [photon/keV]", + 'er_nph': r"$\langle N_{ph}\rangle$ [photon]", + 'er_ne': r"$\langle N_{e^{-}}\rangle$ [electron]", + 'er_nex': r"$\langle N_{ex}\rangle$ [exciton]", + 'er_nion': r"$\langle N_{ion}\rangle$ [ion]", + 'alpha_total': r"$Y_q$ [quanta/keV]", + 'alpha_charge': r"$Y_{e^-}$ [electron/keV]", + 'alpha_light': r"$Y_{\gamma}$ [photon/keV]", + 'alpha_nph': r"$\langle N_{ph}\rangle$ [photon]", + 'alpha_ne': r"$\langle N_{e^{-}}\rangle$ [electron]", + 'alpha_nex': r"$\langle N_{ex}\rangle$ [exciton]", + 'alpha_nion': r"$\langle N_{ion}\rangle$ [ion]", + } + self.titles = { + 'nr_total': "LArNEST Nuclear Recoil Total Yields", + 'nr_charge': "LArNEST Nuclear Recoil Charge Yields", + 'nr_light': "LArNEST Nuclear Recoil Light Yields", + 'nr_nph': "LArNEST Nuclear Recoil Photon Yields", + 'nr_ne': "LArNEST Nuclear Recoil Electron Yields", + 'nr_nex': "LArNEST Nuclear Recoil Exciton Yields", + 'nr_nion': "LArNEST Nuclear Recoil Ion Yields", + 'er_total': "LArNEST Electron Recoil Total Yields", + 'er_charge': "LArNEST Electron Recoil Charge Yields", + 'er_light': "LArNEST Electron Recoil Light Yields", + 'er_nph': "LArNEST Electron Recoil Photon Yields", + 'er_ne': "LArNEST Electron Recoil Electron Yields", + 'er_nex': "LArNEST Electron Recoil Exciton Yields", + 'er_nion': "LArNEST Electron Recoil Ion Yields", + 'alpha_total': "LArNEST Alpha Total Yields", + 'alpha_charge': "LArNEST Alpha Charge Yields", + 'alpha_light': "LArNEST Alpha Light Yields", + 'alpha_nph': "LArNEST Alpha Photon Yields", + 'alpha_ne': "LArNEST Alpha Electron Yields", + 'alpha_nex': "LArNEST Alpha Exciton Yields", + 'alpha_nion': "LArNEST Alpha Ion Yields", + } + self.benchmark_types = { + 'nr_total': ['NR','TotalYields'], + 'nr_charge': ['NR','QuantaYields'], + 'nr_light': ['NR','LightYields'], + 'nr_nph': ['NR','Nph'], + 'nr_ne': ['NR','Ne'], + 'nr_nex': ['NR','Nex'], + 'nr_nion': ['NR','Nion'], + 'er_total': ['ER','TotalYields'], + 'er_charge': ['ER','QuantaYields'], + 'er_light': ['ER','LightYields'], + 'er_nph': ['ER','Nph'], + 'er_ne': ['ER','Ne'], + 'er_nex': ['ER','Nex'], + 'er_nion': ['ER','Nion'], + 'alpha_total': ['Alpha','TotalYields'], + 'alpha_charge': ['Alpha','QuantaYields'], + 'alpha_light': ['Alpha','LightYields'], + 'alpha_nph': ['Alpha','Nph'], + 'alpha_ne': ['Alpha','Ne'], + 'alpha_nex': ['Alpha','Nex'], + 'alpha_nion': ['Alpha','Nion'], + } + + # create plotting directory + if not os.path.isdir("plots/mean_yields/"): + os.makedirs("plots/mean_yields") + if not os.path.isdir("plots/data/"): + os.makedirs("plots/data") + + def plot_data_grid(self, + dataset_type: str='nr_total' + ): + if (len(self.plot_config[dataset_type].keys()) == 1): + return self.plot_data(dataset_type) + fig, axs = generate_plot_grid( + len(self.plot_config[dataset_type].keys()), + figsize=(20,12) + ) + for ii, efield in enumerate(self.plot_config[dataset_type]): + # plot the benchmark values for the given efield + benchmark_mask = (self.benchmarks['type'] == self.benchmark_types[dataset_type][0]) & (self.benchmarks['efield'] == efield) + energy = self.benchmarks['energy'][benchmark_mask] + yields = self.benchmarks[self.benchmark_types[dataset_type][1]][benchmark_mask] + axs.flat[ii].plot(energy, yields, label=f"NEST - {efield} V/cm") + # plot the corresponding data points for this efield + for jj, dataset in enumerate(self.plot_config[dataset_type][efield]): + if (len(self.plot_config[dataset_type][efield][dataset]) < 4): + for kk, field in enumerate(self.plot_config[dataset_type][efield][dataset]): + dataset_mask = (self.datasets[dataset_type]['Dataset'] == dataset) & (self.datasets[dataset_type]['field'] == field) + dataset_energy = self.datasets[dataset_type]['energy'][dataset_mask] + dataset_energy_sl = self.datasets[dataset_type]['energy_sl'][dataset_mask] + dataset_energy_sh = self.datasets[dataset_type]['energy_sh'][dataset_mask] + dataset_yields = self.datasets[dataset_type]['yield'][dataset_mask] + dataset_yields_sl = self.datasets[dataset_type]['yield_sl'][dataset_mask] + dataset_yields_sh = self.datasets[dataset_type]['yield_sh'][dataset_mask] + axs.flat[ii].errorbar( + dataset_energy, dataset_yields, + xerr=[dataset_energy_sl, dataset_energy_sh], + yerr=[dataset_yields_sl, dataset_yields_sh], + linestyle='', marker='.', + label=f"{dataset} - {field} V/cm" + ) + else: + dataset_mask = (self.datasets[dataset_type]['Dataset'] == str(dataset)) & np.isin(np.array(self.datasets[dataset_type]['field']),self.plot_config[dataset_type][efield][dataset]) + dataset_energy = self.datasets[dataset_type]['energy'][dataset_mask] + dataset_energy_sl = self.datasets[dataset_type]['energy_sl'][dataset_mask] + dataset_energy_sh = self.datasets[dataset_type]['energy_sh'][dataset_mask] + dataset_yields = self.datasets[dataset_type]['yield'][dataset_mask] + dataset_yields_sl = self.datasets[dataset_type]['yield_sl'][dataset_mask] + dataset_yields_sh = self.datasets[dataset_type]['yield_sh'][dataset_mask] + axs.flat[ii].errorbar( + dataset_energy, dataset_yields, + xerr=[dataset_energy_sl, dataset_energy_sh], + yerr=[dataset_yields_sl, dataset_yields_sh], + linestyle='', marker='.', + label=f"{dataset} - {min(self.plot_config[dataset_type][efield][dataset])}-{max(self.plot_config[dataset_type][efield][dataset])} V/cm" + ) + axs.flat[ii].set_xscale("log") + if ("_n" in dataset_type): + axs.flat[ii].set_yscale("log") + axs.flat[ii].legend() + fig.supxlabel("Energy [keV]") + fig.supylabel(self.ylabels[dataset_type]) + plt.suptitle(self.titles[dataset_type]) + plt.tight_layout() + plt.savefig(f"plots/data/{self.benchmark_types[dataset_type][0]}_{self.benchmark_types[dataset_type][1]}_Data.png") + plt.close() + + def plot_data(self, + dataset_type: str='nr_total' + ): + for ii, efield in enumerate(self.plot_config[dataset_type]): + fig, axs = plt.subplots() + # plot the benchmark values for the given efield + benchmark_mask = (self.benchmarks['type'] == self.benchmark_types[dataset_type][0]) & (self.benchmarks['efield'] == efield) + energy = self.benchmarks['energy'][benchmark_mask] + yields = self.benchmarks[self.benchmark_types[dataset_type][1]][benchmark_mask] + axs.plot(energy, yields, label=f"NEST - {efield} V/cm") + # plot the corresponding data points for this efield + for jj, dataset in enumerate(self.plot_config[dataset_type][efield]): + if (len(self.plot_config[dataset_type][efield][dataset]) < 4): + for kk, field in enumerate(self.plot_config[dataset_type][efield][dataset]): + dataset_mask = (self.datasets[dataset_type]['Dataset'] == str(dataset)) & (self.datasets[dataset_type]['field'] == field) + dataset_energy = self.datasets[dataset_type]['energy'][dataset_mask] + dataset_energy_sl = self.datasets[dataset_type]['energy_sl'][dataset_mask] + dataset_energy_sh = self.datasets[dataset_type]['energy_sh'][dataset_mask] + dataset_yields = self.datasets[dataset_type]['yield'][dataset_mask] + dataset_yields_sl = self.datasets[dataset_type]['yield_sl'][dataset_mask] + dataset_yields_sh = self.datasets[dataset_type]['yield_sh'][dataset_mask] + axs.errorbar( + dataset_energy, dataset_yields, + xerr=[dataset_energy_sl, dataset_energy_sh], + yerr=[dataset_yields_sl, dataset_yields_sh], + linestyle='', marker='.', + label=f"{dataset} - {field} V/cm" + ) + else: + dataset_mask = (self.datasets[dataset_type]['Dataset'] == str(dataset)) & np.isin(np.array(self.datasets[dataset_type]['field']),self.plot_config[dataset_type][efield][dataset]) + dataset_energy = self.datasets[dataset_type]['energy'][dataset_mask] + dataset_energy_sl = self.datasets[dataset_type]['energy_sl'][dataset_mask] + dataset_energy_sh = self.datasets[dataset_type]['energy_sh'][dataset_mask] + dataset_yields = self.datasets[dataset_type]['yield'][dataset_mask] + dataset_yields_sl = self.datasets[dataset_type]['yield_sl'][dataset_mask] + dataset_yields_sh = self.datasets[dataset_type]['yield_sh'][dataset_mask] + axs.errorbar( + dataset_energy, dataset_yields, + xerr=[dataset_energy_sl, dataset_energy_sh], + yerr=[dataset_yields_sl, dataset_yields_sh], + linestyle='', marker='.', + label=f"{dataset} - {min(self.plot_config[dataset_type][efield][dataset])}-{max(self.plot_config[dataset_type][efield][dataset])} V/cm" + ) + axs.set_xscale("log") + if ("_n" in dataset_type): + axs.set_yscale("log") + plt.legend(bbox_to_anchor=(1.05, 1.0), loc='upper left') + axs.set_xlabel("Energy [keV]") + axs.set_ylabel(self.ylabels[dataset_type]) + plt.title(self.titles[dataset_type]) + plt.tight_layout() + plt.savefig(f"plots/data/{self.benchmark_types[dataset_type][0]}_{self.benchmark_types[dataset_type][1]}_{efield}_Data.png") + plt.close() + + def plot_yields(self, + dataset_type: str='nr_total' + ): + fig, axs = plt.subplots(figsize=(10,6)) + # plot the benchmark values for the given efield + benchmark_mask = (self.benchmarks['type'] == self.benchmark_types[dataset_type][0]) + unique_efields = np.unique(self.benchmarks['efield'][benchmark_mask]) + for efield in unique_efields: + efield_mask = (self.benchmarks['efield'][benchmark_mask] == efield) + energy = self.benchmarks['energy'][benchmark_mask][efield_mask] + yields = self.benchmarks[self.benchmark_types[dataset_type][1]][benchmark_mask][efield_mask] + axs.plot(energy, yields, label=f"NEST - {efield} V/cm") + axs.set_xscale("log") + if ("_n" in dataset_type): + axs.set_yscale("log") + plt.legend(bbox_to_anchor=(1.05, 1.0), loc='upper left') + axs.set_xlabel("Energy [keV]") + axs.set_ylabel(self.ylabels[dataset_type]) + plt.title(self.titles[dataset_type]) + plt.tight_layout() + plt.savefig(f"plots/mean_yields/{self.benchmark_types[dataset_type][0]}_{self.benchmark_types[dataset_type][1]}.png") + plt.close() + + def plot_all_yields(self, + ): + for benchmark_type in self.benchmark_types: + self.plot_yields(benchmark_type) \ No newline at end of file diff --git a/utils/larnest/parameters.py b/utils/larnest/parameters.py new file mode 100644 index 00000000..71d5284a --- /dev/null +++ b/utils/larnest/parameters.py @@ -0,0 +1,166 @@ +""" +Collections of datasets and information. +""" +nr_total_yield_default = { + 100: { + "SCENE 2015": [96.4, 193, 293], + "ARIS 2018": [200], + }, +} +nr_charge_yield_default = { + 100: { + "SCENE 2015": [96.4], + }, + 250: { + "SCENE 2015": [193, 293], + "ARIS 2018": [200], + "Joshi 2014": [240], + }, + 500: { + "SCENE 2015": [486], + "Joshi 2014": [640], + }, + 1750: { + "Joshi 2014": [1600, 2130], + }, +} +nr_light_yield_default = { + 1: { + "ARIS 2018": [0], + "SCENE 2015": [0], + "Regenfus 2012":[0], + "WARP 2005": [0], + "CREUS 2015": [0], + "MicroClean 2012": [0], + }, + 50: { + "SCENE 2013": [50], + "ARIS 2018": [50], + }, + 100: { + "SCENE 2015": [100], + "ARIS 2018": [100], + }, + 250: { + "SCENE 2015": [193, 293], + "ARIS 2018": [200], + "SCENE 2013": [300], + }, + 500: { + "ARIS 2018": [500], + }, +} +er_charge_yield_default = { + 1: { + "ARIS 2018": [0], + "DarkSide10 2013": [0], + "Lippincott 2010": [0], + "WARP 2005": [0], + "Kimura 2019": [0], + "Doke 2002": [0] + }, + 100: { + "ARIS 2018": [50,100], + "Scalettar 1982": [84,89,94,101,110,119,128,129,130,139,145,148], + }, + 200: { + "ARIS 2018": [200], + "Joshi 2014": [200], + "Scalettar 1982":[152,160,176,180,185,188,192,200,205,210,240,243,248,268,270,276,285,309,313,351,371,375,388,391] + }, + 600: { + "ARIS 2018": [500], + "Joshi 2014": [550], + "Scalettar 1982": [410,411,427,441,478,481,490,510,531,536,601,620,626,661], + "Bondar 2016": [600], + #"Aprile ": [572] + }, + 1000: { + "Joshi 2014": [1200], + "Scalettar 1982": [801,820,832,899,904,943], + #"Aprile ": [] + }, + 1500: { + "Joshi 2014": [1600,1750], + "Scalettar 1982": [1004,1012,1064,1244,1403,1455], + "Bondar 2016": [1750], + }, + 2500: { + "Joshi 2014": [2150,2400,3000], + "Scalettar 1982":[2009,2064,2806,2913], + "Bondar 2016": [2400], + "Sangiorgio": [2400], + "Doke 2002": [2010], + }, + 6000: { + "Scalettar 1982":[4600,5704,6586,6693], + "Doke 2002": [4020,5000,6010], + #"Aprile", + }, + 9500: { + "Scalettar 1982":[8490,8673,9481,9681], + "Doke 2002": [8000,9000,9990], + #"Aprile", + }, +} +er_light_yield_default = { + 1: { + "ARIS 2018": [0], + "DarkSide10 2013": [0], + "Lippincott 2010": [0], + "WARP 2005": [0], + "Kimura 2019": [0], + "Kimura LIDINE":[0], + }, + 100: { + "ARIS 2018": [50,100], + "Scalettar 1982": [84,89,94,101,110,119,128,129,130,139,145,148], + }, + 200: { + "ARIS 2018": [200], + "Joshi 2014": [200], + "Scalettar 1982":[152,160,176,180,185,188,192,200,205,210,240,243,248,268,270,276,285,309,313,351,371,375,388,391] + }, + 600: { + "ARIS 2018": [500], + "Joshi 2014": [550], + "Scalettar 1982": [410,411,427,441,478,481,490,510,531,536,601,620,626,661], + "Bondar 2016": [600], + #"Aprile ": [572] + }, + 1000: { + "Joshi 2014": [1200], + "Scalettar 1982": [801,820,832,899,904,943], + #"Aprile ": [] + }, + 1500: { + "Joshi 2014": [1600,1750], + "Scalettar 1982": [1004,1012,1064,1244,1403,1455], + "Bondar 2016": [1750], + }, + 2500: { + "Joshi 2014": [2150,2400,3000], + "Scalettar 1982":[2009,2064,2806,2913], + "Bondar 2016": [2400], + "Sangiorgio": [2400], + "Doke 2002": [2010], + }, + 6000: { + "Scalettar 1982":[4600,5704,6586,6693], + "Doke 2002": [4020,5000,6010], + #"Aprile", + }, + 9500: { + "Scalettar 1982":[8490,8673,9481,9681], + "Doke 2002": [8000,9000,9990], + #"Aprile", + }, +} + +default_plot_config = { + "nr_total": nr_total_yield_default, + "nr_charge": nr_charge_yield_default, + "nr_light": nr_light_yield_default, + "er_charge": er_charge_yield_default, + "er_light": er_light_yield_default, +} \ No newline at end of file diff --git a/utils/larnest/plot_benchmarks.py b/utils/larnest/plot_benchmarks.py new file mode 100644 index 00000000..95e45dd7 --- /dev/null +++ b/utils/larnest/plot_benchmarks.py @@ -0,0 +1,40 @@ +""" +Benchmark functions for LArNEST +""" +import numpy as np +from matplotlib import pyplot as plt +import pandas as pd +import os + +from parameters import * +from lar_dataset import LArDataset + +if __name__ == "__main__": + + dataset_file = "data/lar_data.npz" + benchmarks_file = "data/benchmarks.csv" + + # create the dataset object + lar_dataset = LArDataset( + dataset_file=dataset_file, + benchmarks_file=benchmarks_file, + plot_config=default_plot_config, + ) + + # make benchmark plots + lar_dataset.plot_all_yields() + + # plot with data + # nr yields + lar_dataset.plot_data('nr_total') + lar_dataset.plot_data_grid('nr_total') + lar_dataset.plot_data('nr_charge') + lar_dataset.plot_data_grid('nr_charge') + lar_dataset.plot_data('nr_light') + lar_dataset.plot_data_grid('nr_light') + + # er yields + lar_dataset.plot_data('er_charge') + lar_dataset.plot_data_grid('er_charge') + lar_dataset.plot_data('er_light') + lar_dataset.plot_data_grid('er_light') \ No newline at end of file diff --git a/utils/larnest/utils.py b/utils/larnest/utils.py new file mode 100644 index 00000000..8141533d --- /dev/null +++ b/utils/larnest/utils.py @@ -0,0 +1,174 @@ +""" +Various utilities for plotting, etc. +""" +import numpy as np +from matplotlib import pyplot as plt +import csv + + +def generate_plot_grid( + num_plots, + **kwargs, +): + """ + Generate a grid of N plots, excluding extra + ones which are not used. + """ + nrows = int(np.floor(np.sqrt(num_plots))) + ncols = int(np.ceil(num_plots/nrows)) + fig, axs = plt.subplots( + nrows=nrows, ncols=ncols, + **kwargs + ) + nplots = nrows * ncols + nextra = nplots - num_plots + for ii in range(nextra): + axs.flat[-(ii+1)].set_visible(False) + return fig, axs + +if __name__ == "__main__": + + nr_total = { + "Dataset": [], + "energy": [], + "energy_sl":[], + "energy_sh":[], + "field": [], + "field_sl": [], + "field_sh": [], + "yield": [], + "yield_sl": [], + "yield_sh": [] + } + with open("data/nr_total_alt.csv","r") as file: + reader = csv.reader(file, delimiter=",") + next(reader) + for row in reader: + nr_total["Dataset"].append(row[0]) + nr_total["energy"].append(float(row[1])) + nr_total["energy_sl"].append(float(row[2])) + nr_total["energy_sh"].append(float(row[3])) + nr_total["field"].append(float(row[4])) + nr_total["field_sl"].append(float(row[5])) + nr_total["field_sh"].append(float(row[6])) + nr_total["yield"].append(float(row[7])) + nr_total["yield_sl"].append(float(row[8])) + nr_total["yield_sh"].append(float(row[9])) + nr_total = {key: np.array(nr_total[key]) for key in nr_total} + nr_charge = { + "Dataset": [], + "energy": [], + "energy_sl":[], + "energy_sh":[], + "field": [], + "field_sl": [], + "field_sh": [], + "yield": [], + "yield_sl": [], + "yield_sh": [] + } + with open("data/nr_charge_alt.csv","r") as file: + reader = csv.reader(file, delimiter=",") + next(reader) + for row in reader: + nr_charge["Dataset"].append(row[0]) + nr_charge["energy"].append(float(row[1])) + nr_charge["energy_sl"].append(float(row[2])) + nr_charge["energy_sh"].append(float(row[3])) + nr_charge["field"].append(float(row[4])) + nr_charge["field_sl"].append(float(row[5])) + nr_charge["field_sh"].append(float(row[6])) + nr_charge["yield"].append(float(row[7])) + nr_charge["yield_sl"].append(float(row[8])) + nr_charge["yield_sh"].append(float(row[9])) + nr_charge = {key: np.array(nr_charge[key]) for key in nr_charge} + nr_light = { + "Dataset": [], + "energy": [], + "energy_sl":[], + "energy_sh":[], + "field": [], + "field_sl": [], + "field_sh": [], + "yield": [], + "yield_sl": [], + "yield_sh": [] + } + with open("data/nr_light_alt.csv","r") as file: + reader = csv.reader(file, delimiter=",") + next(reader) + for row in reader: + nr_light["Dataset"].append(row[0]) + nr_light["energy"].append(float(row[1])) + nr_light["energy_sl"].append(float(row[2])) + nr_light["energy_sh"].append(float(row[3])) + nr_light["field"].append(float(row[4])) + nr_light["field_sl"].append(float(row[5])) + nr_light["field_sh"].append(float(row[6])) + nr_light["yield"].append(float(row[7])) + nr_light["yield_sl"].append(float(row[8])) + nr_light["yield_sh"].append(float(row[9])) + nr_light = {key: np.array(nr_light[key]) for key in nr_light} + er_charge = { + "Dataset": [], + "energy": [], + "energy_sl":[], + "energy_sh":[], + "field": [], + "field_sl": [], + "field_sh": [], + "yield": [], + "yield_sl": [], + "yield_sh": [] + } + with open("data/er_charge.csv","r") as file: + reader = csv.reader(file, delimiter=",") + next(reader) + for row in reader: + er_charge["Dataset"].append(row[0]) + er_charge["energy"].append(float(row[1])) + er_charge["energy_sl"].append(float(row[2])) + er_charge["energy_sh"].append(float(row[3])) + er_charge["field"].append(float(row[4])) + er_charge["field_sl"].append(float(row[5])) + er_charge["field_sh"].append(float(row[6])) + er_charge["yield"].append(float(row[7])) + er_charge["yield_sl"].append(float(row[8])) + er_charge["yield_sh"].append(float(row[9])) + er_charge = {key: np.array(er_charge[key]) for key in er_charge} + er_light = { + "Dataset": [], + "energy": [], + "energy_sl":[], + "energy_sh":[], + "field": [], + "field_sl": [], + "field_sh": [], + "yield": [], + "yield_sl": [], + "yield_sh": [] + } + with open("data/er_light.csv","r") as file: + reader = csv.reader(file, delimiter=",") + next(reader) + for row in reader: + er_light["Dataset"].append(row[0]) + er_light["energy"].append(float(row[1])) + er_light["energy_sl"].append(float(row[2])) + er_light["energy_sh"].append(float(row[3])) + er_light["field"].append(float(row[4])) + er_light["field_sl"].append(float(row[5])) + er_light["field_sh"].append(float(row[6])) + er_light["yield"].append(float(row[7])) + er_light["yield_sl"].append(float(row[8])) + er_light["yield_sh"].append(float(row[9])) + er_light = {key: np.array(er_light[key]) for key in er_light} + + np.savez( + "data/lar_data.npz", + nr_total=nr_total, + nr_charge=nr_charge, + nr_light=nr_light, + er_charge=er_charge, + er_light=er_light, + ) \ No newline at end of file