diff --git a/examples/LArNEST/CMakeLists.txt b/examples/LArNEST/CMakeLists.txt index 2db85ae6..c433ee04 100644 --- a/examples/LArNEST/CMakeLists.txt +++ b/examples/LArNEST/CMakeLists.txt @@ -4,5 +4,11 @@ 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 +add_executable(LArNESTMeanYieldsBenchmarks LArNESTMeanYieldsBenchmarks.cpp) +target_link_libraries(LArNESTMeanYieldsBenchmarks PUBLIC NEST::Core) + +add_executable(LArNESTFluctuationBenchmarks LArNESTFluctuationBenchmarks.cpp) +target_link_libraries(LArNESTFluctuationBenchmarks PUBLIC NEST::Core) + +add_executable(LArNESTNeutronCapture LArNESTNeutronCapture.cpp) +target_link_libraries(LArNESTNeutronCapture PUBLIC NEST::Core) \ No newline at end of file diff --git a/examples/LArNEST/LArNESTFluctuationBenchmarks.cpp b/examples/LArNEST/LArNESTFluctuationBenchmarks.cpp new file mode 100644 index 00000000..1c2709c8 --- /dev/null +++ b/examples/LArNEST/LArNESTFluctuationBenchmarks.cpp @@ -0,0 +1,173 @@ +/** + * @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; + uint64_t num_events = 100; + if (argc > 1) { + density = atof(argv[1]); + } + if (argc > 2) { + seed = atoi(argv[2]); + } + if (argc > 3) { + num_events = atoi(argv[3]); + } + + // 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("fluctuation_benchmarks.csv"); + output_file << "type,energy,efield,TotalYield,QuantaYield,LightYield,Nph,Ne,Nex,Nion,TotalYield_std,QuantaYield_std,LightYield_std,Nph_std,Ne_std,Nex_std,Nion_std\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 + }); + } + // 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 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.FullCalculation( + particle_types[k], + energy_vals[i], + electric_field[v], + density, + false + ); + Nph[j] = double(result.fluctuations.NphFluctuation); + Ne[j] = double(result.fluctuations.NeFluctuation); + Nex[j] = double(result.fluctuations.NexFluctuation); + Nion[j] = double(result.fluctuations.NionFluctuation); + } + 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); + + double Nq_mean = Nph_mean + Ne_mean; + double Nq_std = 0.0; + double Nph_std = 0.0; + double Ne_std = 0.0; + double Nex_std = 0.0; + double Nion_std = 0.0; + + for (size_t ii = 0; ii < Nph.size(); ii++) { + Nq_std += (Nph[ii] + Ne[ii] - Nq_mean) * (Nph[ii] + Ne[ii] - Nq_mean); + } + Nq_std = std::sqrt(Nq_std / double(num_events)); + + for (size_t ii = 0; ii < Nph.size(); ii++) { + Nph_std += (Nph[ii] - Nph_mean) * (Nph[ii] - Nph_mean); + } + Nph_std = std::sqrt(Nph_std / double(num_events)); + + for (size_t ii = 0; ii < Ne.size(); ii++) { + Ne_std += (Ne[ii] - Ne_mean) * (Ne[ii] - Ne_mean); + } + Ne_std = std::sqrt(Ne_std / double(num_events)); + + for (size_t ii = 0; ii < Nex.size(); ii++) { + Nex_std += (Nex[ii] - Nex_mean) * (Nex[ii] - Nex_mean); + } + Nex_std = std::sqrt(Nex_std / double(num_events)); + + for (size_t ii = 0; ii < Nion.size(); ii++) { + Nion_std += (Nion[ii] - Nion_mean) * (Nion[ii] - Nion_mean); + } + Nion_std = std::sqrt(Nion_std / double(num_events)); + + output_file << particle_type[k] << ","; + output_file << energy_vals[i] << ","; + output_file << electric_field[v] << ","; + output_file << (Nq_mean / energy_vals[i]) << ","; + output_file << (Ne_mean / energy_vals[i]) << ","; + output_file << (Nph_mean / energy_vals[i]) << ","; + output_file << Nph_mean << ","; + output_file << Ne_mean << ","; + output_file << Nex_mean << ","; + output_file << Nion_mean << ","; + output_file << (Nq_std / energy_vals[i]) << ","; + output_file << (Ne_std / energy_vals[i]) << ","; + output_file << (Nph_std / energy_vals[i]) << ","; + output_file << Nph_std << ","; + output_file << Ne_std << ","; + output_file << Nex_std << ","; + output_file << Nion_std << "\n"; + } + } + } + output_file.close(); + return 0; +} \ No newline at end of file diff --git a/examples/LArNEST/LArNESTBenchmarks.cpp b/examples/LArNEST/LArNESTMeanYieldsBenchmarks.cpp similarity index 93% rename from examples/LArNEST/LArNESTBenchmarks.cpp rename to examples/LArNEST/LArNESTMeanYieldsBenchmarks.cpp index 1d26e562..60eeed6a 100644 --- a/examples/LArNEST/LArNESTBenchmarks.cpp +++ b/examples/LArNEST/LArNESTMeanYieldsBenchmarks.cpp @@ -58,8 +58,8 @@ int main(int argc, char* argv[]) 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"; + output_file.open("mean_yields_benchmarks.csv"); + output_file << "type,energy,efield,TotalYield,QuantaYield,LightYield,Nph,Ne,Nex,Nion,TotalYield_std,QuantaYield_std,LightYield_std,Nph_std,Ne_std,Nex_std,Nion_std\n"; // iterate over electric field values for (size_t k = 0; k < particle_types.size(); k++) @@ -104,7 +104,7 @@ int main(int argc, char* argv[]) output_file << result.yields.Nph << ","; output_file << result.yields.Ne << ","; output_file << result.yields.Nex << ","; - output_file << result.yields.Nion << "\n"; + output_file << result.yields.Nion << ",0,0,0,0,0,0,0\n"; } } } diff --git a/examples/LArNEST/LArNESTNeutronCapture.cpp b/examples/LArNEST/LArNESTNeutronCapture.cpp new file mode 100644 index 00000000..62b16db5 --- /dev/null +++ b/examples/LArNEST/LArNESTNeutronCapture.cpp @@ -0,0 +1,84 @@ +/** + * @file LArNESTNeutronCapture.cpp + * @author Nicholas Carrara [nmcarrara@ucdavis.edu] + * @brief + * @version 0.1 + * @date 2022-05-23 + */ +#include +#include +#include +#include +#include + +#include "LArDetector.hh" +#include "LArNEST.hh" + +int main(int argc, char* argv[]) +{ + double density = 1.393; + uint64_t seed = 0; + if (argc > 1) { + density = atof(argv[1]); + } + if (argc > 2) { + seed = atoi(argv[2]); + } + + std::vector electric_field = { + 1, 100, 200, 500, 600, 1000, 1500, 2500, 6000, 9000, 9500 + }; + + std::vector energy_vals = { + 167.3, 516.0, 348.7, 867.3, 837.7, 1186.8, 1354.0, + 1044.3, 1881.5, 2229.5, 2566.1, 2432.5, 2781.8, + 2842.5, 3111.3, 1972.6, 2291.6, 2810.5, + 3564.7, 3405.5, 2668.1, 3564.5, 2614.3, + 3451.8, 4102.5, + 1828.8, 2130.7, 2668.1, 2771.8, 3089.4, 3150.2, + 3365.5, 3405.3, 3700.4, 4745.0, 5063.7, 5582.0, + 6098.9 + }; + + 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("neutron_capture.csv"); + output_file << "type,energy,efield,TotalYield,QuantaYield,LightYield,Nph,Ne,Nex,Nion\n"; + + for (size_t v = 0; v < electric_field.size(); v++) + { + // iterate over energy values + for (size_t i = 0; i < energy_vals.size(); i++) + { + result = larnest.FullCalculation( + NEST::LArInteraction::ER, + energy_vals[i], + electric_field[v], + density, + false + ); + output_file << "ER,"; + 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 index cbc8abf0..38beb008 100644 --- a/examples/LArNEST/LegacyLArNESTBenchmarks.cpp +++ b/examples/LArNEST/LegacyLArNESTBenchmarks.cpp @@ -63,7 +63,7 @@ int main(int argc, char* argv[]) // 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.open("legacy_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++) diff --git a/include/NEST/LArNEST.hh b/include/NEST/LArNEST.hh index 32c17225..3b4592ae 100644 --- a/include/NEST/LArNEST.hh +++ b/include/NEST/LArNEST.hh @@ -32,6 +32,11 @@ namespace NEST Alpha = 2, }; + enum class LArFluctuationModel + { + Default = 0, + }; + struct LArNRYieldsParameters { double alpha = {11.10}; @@ -175,9 +180,6 @@ namespace NEST struct LArYieldFluctuationResult { - double TotalYieldFluctuation; - double QuantaYieldFluctuation; - double LightYieldFluctuation; double NphFluctuation; double NeFluctuation; double NexFluctuation; @@ -238,15 +240,15 @@ namespace NEST 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; } + 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; } @@ -275,7 +277,8 @@ namespace NEST * */ LArYieldFluctuationResult GetYieldFluctuations( - const YieldResult &yields, double density + const LArYieldResult &yields, + double density ); /** * @brief @@ -443,7 +446,17 @@ namespace NEST LArYieldResult GetAlphaYields( double energy, double efield, double density ); - + //-------------------------Fluctuation Yields-------------------------// + /** + * @brief Get the Default Fluctuations object + * + * @param yields + * @param density + * @return LArYieldFluctuationResult + */ + LArYieldFluctuationResult GetDefaultFluctuations( + const LArYieldResult &yields, double density + ); //-------------------------Photon Times-------------------------// double GetPhotonTime( LArInteraction species, bool exciton, @@ -518,6 +531,7 @@ namespace NEST double fWorkPhotonFunction = {14.544}; double fNexOverNion = {0.21}; + double fALF = {1. / (1. + 0.21)}; double fFanoER = {0.1115}; LArNRYieldsParameters fNR; @@ -526,5 +540,7 @@ namespace NEST ThomasImelParameters fThomasImelParameters; DriftParameters fDriftParameters; + + enum LArFluctuationModel fLArFluctuationModel; }; } diff --git a/src/LArNEST.cpp b/src/LArNEST.cpp index 299b4168..ccd0c821 100644 --- a/src/LArNEST.cpp +++ b/src/LArNEST.cpp @@ -27,14 +27,12 @@ namespace NEST 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; - } + // Thomas-Imel is not used here, since alpha = Nex/Nion + // is assumed to be a fixed quantity, and so the + // recombination probability is constrained through + // the yields functions: + // Nion = Nq/(1 + alpha) = (Ye + Yph)E/(1 + alpha) + double Nion = Nq / (1 + fNexOverNion); double Nex = Nq - Nion; double WQ_eV = fWorkQuantaFunction; @@ -62,12 +60,16 @@ namespace NEST } LArYieldFluctuationResult LArNEST::GetYieldFluctuations( - const YieldResult &yields, double density + const LArYieldResult &yields, + double density ) { - // TODO: Understand and break up this function. - LArYieldFluctuationResult result{}; - return result; // quanta returned with recomb fluctuations + if (fLArFluctuationModel == LArFluctuationModel::Default) { + return GetDefaultFluctuations(yields, density); + } + else { + return GetDefaultFluctuations(yields, density); + } } LArNESTResult LArNEST::FullCalculation( @@ -79,7 +81,7 @@ namespace NEST if (density < 1.) fdetector->set_inGas(true); LArNESTResult result; result.yields = GetYields(species, energy, efield, density); - // result.quanta = GetQuanta(result.yields, density); + result.fluctuations = GetYieldFluctuations(result.yields, density); // if (do_times) // result.photon_times = GetPhotonTimes( // species, result.quanta.photons, @@ -247,6 +249,69 @@ namespace NEST energy, efield ); } + //-------------------------Fluctuation Yields-------------------------// + LArYieldFluctuationResult LArNEST::GetDefaultFluctuations( + const LArYieldResult &yields, double density + ) + { + double Fano = 0.1115; + double NexOverNion = yields.Nex/yields.Nion; + double p_r = (yields.Nph - yields.Nex)/yields.Nion; + double alf = 1./(1. + NexOverNion); + double Nq = 0; + double Nex = 0; + double Nion = 0; + double Nq_mean = yields.Ne + yields.Nph; + if (yields.Lindhard == 1.) + { + Nq = floor( + RandomGen::rndm()->rand_gauss( + Nq_mean, + sqrt(Fano * Nq_mean), + true + ) + 0.5 + ); + if (Nq < 0.) { + Nq = 0.; + } + Nion = double(RandomGen::rndm()->binom_draw(Nq, alf)); + if (Nion > Nq) { + Nion = Nq; + } + Nex = Nq - Nion; + } + else + { + Nion = floor( + RandomGen::rndm()->rand_gauss( + Nq_mean * alf, + sqrt(Fano * Nq_mean * alf), + true + ) + 0.5 + ); + Nex = floor( + RandomGen::rndm()->rand_gauss( + Nq_mean * alf * NexOverNion, + sqrt(Fano * Nq_mean * alf * NexOverNion), + true + ) + 0.5 + ); + if (Nex < 0.) { + Nex = 0.; + } + if (Nion < 0.) { + Nion = 0.; + } + Nq = Nion + Nex; + } + LArYieldFluctuationResult result{ + Nex + p_r*Nion, + (1. - p_r)*Nion, + Nex, + Nion + }; + return result; + } //-------------------------Photon Times-------------------------// double LArNEST::GetPhotonTime( @@ -471,6 +536,7 @@ namespace NEST 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) { @@ -534,7 +600,7 @@ namespace NEST ) { // this section calculates recombination following the modified - // Birks'Law of Doke, deposition by deposition, may be overridden + // 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; diff --git a/utils/larnest/lar_dataset.py b/utils/larnest/lar_dataset.py index b5cf19cd..6a0d2555 100644 --- a/utils/larnest/lar_dataset.py +++ b/utils/larnest/lar_dataset.py @@ -15,13 +15,15 @@ class LArDataset: """ """ def __init__(self, - dataset_file: str='data/lar_data.npz', + dataset_file: str='data/lar_data.npz', benchmarks_file: str='', - plot_config: dict=default_plot_config, + plot_config: dict=default_plot_config, + fluctuations: bool=False, ): self.dataset_file = dataset_file self.benchmarks_file = benchmarks_file self.plot_config = plot_config + self.fluctuations = fluctuations # 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} @@ -31,7 +33,9 @@ def __init__(self, names=[ 'type','energy', 'efield' , 'TotalYields', 'QuantaYields', 'LightYields', - 'Nph' , 'Ne', 'Nex', 'Nion' + 'Nph' , 'Ne', 'Nex', 'Nion', + 'TotalYields_std', 'QuantaYields_std', 'LightYields_std', + 'Nph_std', 'Ne_std', 'Nex_std', 'Nion_std', ], header=0, #dtype=float @@ -107,11 +111,19 @@ def __init__(self, 'alpha_nion': ['Alpha','Nion'], } + self.fluctuation_types = [ + 'Nph', 'Ne', 'Nex', 'Nion' + ] + # create plotting directory if not os.path.isdir("plots/mean_yields/"): os.makedirs("plots/mean_yields") + if not os.path.isdir("plots/fluctuations/"): + os.makedirs("plots/fluctuations") if not os.path.isdir("plots/data/"): os.makedirs("plots/data") + if not os.path.isdir("plots/data_fluctuations/"): + os.makedirs("plots/data_fluctuations") def plot_data_grid(self, dataset_type: str='nr_total' @@ -127,7 +139,14 @@ def plot_data_grid(self, 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] + yields[(yields<0)] = 0 axs.flat[ii].plot(energy, yields, label=f"NEST - {efield} V/cm") + if self.fluctuations: + errors = self.benchmarks[self.benchmark_types[dataset_type][1]+"_std"][benchmark_mask] + low_errors = yields - errors + high_errors = yields + errors + low_errors[(low_errors)<0] = 0 + axs.flat[ii].fill_between(energy, low_errors, high_errors) # 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): @@ -169,7 +188,10 @@ def plot_data_grid(self, 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") + if self.fluctuations: + plt.savefig(f"plots/data_fluctuations/{self.benchmark_types[dataset_type][0]}_{self.benchmark_types[dataset_type][1]}_Data.png") + else: + 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, @@ -181,7 +203,14 @@ def plot_data(self, 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] + yields[(yields<0)] = 0 axs.plot(energy, yields, label=f"NEST - {efield} V/cm") + if self.fluctuations: + errors = self.benchmarks[self.benchmark_types[dataset_type][1]+"_std"][benchmark_mask] + low_errors = yields - errors + high_errors = yields + errors + low_errors[(low_errors)<0] = 0 + axs.fill_between(energy, low_errors, high_errors) # 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): @@ -223,10 +252,13 @@ def plot_data(self, 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") + if self.fluctuations: + plt.savefig(f"plots/data_fluctuations/{self.benchmark_types[dataset_type][0]}_{self.benchmark_types[dataset_type][1]}_{efield}_Data.png") + else: + 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, + def plot_mean_yields(self, dataset_type: str='nr_total' ): fig, axs = plt.subplots(figsize=(10,6)) @@ -237,6 +269,7 @@ def plot_yields(self, 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] + yields[(yields<0)] = 0 axs.plot(energy, yields, label=f"NEST - {efield} V/cm") axs.set_xscale("log") if ("_n" in dataset_type): @@ -252,4 +285,38 @@ def plot_yields(self, def plot_all_yields(self, ): for benchmark_type in self.benchmark_types: - self.plot_yields(benchmark_type) \ No newline at end of file + self.plot_mean_yields(benchmark_type) + + def plot_fluctuations(self, + dataset_type: str='nr_total' + ): + # 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: + fig, axs = plt.subplots(figsize=(10,6)) + 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] + yields[(yields<0)] = 0 + errors = self.benchmarks[self.benchmark_types[dataset_type][1]+"_std"][benchmark_mask][efield_mask] + low_errors = yields - errors + high_errors = yields + errors + low_errors[(low_errors)<0] = 0 + axs.fill_between(energy, low_errors, high_errors, color='y', alpha=0.7) + axs.plot(energy, yields, label=f"NEST - {efield} V/cm", color='k', linestyle='-') + 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/fluctuations/{self.benchmark_types[dataset_type][0]}_{self.benchmark_types[dataset_type][1]}_{efield}Vcm.png") + plt.close() + + def plot_all_fluctuations(self, + ): + for benchmark_type in self.benchmark_types: + self.plot_fluctuations(benchmark_type) \ No newline at end of file diff --git a/utils/larnest/plot_benchmarks.py b/utils/larnest/plot_benchmarks.py index 95e45dd7..64945be7 100644 --- a/utils/larnest/plot_benchmarks.py +++ b/utils/larnest/plot_benchmarks.py @@ -12,29 +12,40 @@ if __name__ == "__main__": dataset_file = "data/lar_data.npz" - benchmarks_file = "data/benchmarks.csv" - + mean_yields_benchmarks_file = "data/mean_yields_benchmarks.csv" + fluctuations_benchmarks_file = "data/fluctuation_benchmarks.csv" + lar_datasets = [] + # create the dataset object - lar_dataset = LArDataset( + lar_datasets.append(LArDataset( + dataset_file=dataset_file, + benchmarks_file=mean_yields_benchmarks_file, + plot_config=default_plot_config, + fluctuations=False, + )) + lar_datasets.append(LArDataset( dataset_file=dataset_file, - benchmarks_file=benchmarks_file, + benchmarks_file=fluctuations_benchmarks_file, plot_config=default_plot_config, - ) + fluctuations=True, + )) # make benchmark plots - lar_dataset.plot_all_yields() + lar_datasets[0].plot_all_yields() + lar_datasets[1].plot_all_fluctuations() - # 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') + for ii, lar_dataset in enumerate(lar_datasets): + # 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 + # 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 index 8141533d..ef3a1f34 100644 --- a/utils/larnest/utils.py +++ b/utils/larnest/utils.py @@ -26,8 +26,9 @@ def generate_plot_grid( axs.flat[-(ii+1)].set_visible(False) return fig, axs -if __name__ == "__main__": - +def compile_nr_total_data( + data_dir: str="data/", +): nr_total = { "Dataset": [], "energy": [], @@ -40,7 +41,7 @@ def generate_plot_grid( "yield_sl": [], "yield_sh": [] } - with open("data/nr_total_alt.csv","r") as file: + with open(f"{data_dir}nr_total_alt.csv","r") as file: reader = csv.reader(file, delimiter=",") next(reader) for row in reader: @@ -54,7 +55,11 @@ def generate_plot_grid( 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} + return {key: np.array(nr_total[key]) for key in nr_total} + +def compile_nr_charge_data( + data_dir: str="data/", +): nr_charge = { "Dataset": [], "energy": [], @@ -67,7 +72,7 @@ def generate_plot_grid( "yield_sl": [], "yield_sh": [] } - with open("data/nr_charge_alt.csv","r") as file: + with open(f"{data_dir}nr_charge_alt.csv","r") as file: reader = csv.reader(file, delimiter=",") next(reader) for row in reader: @@ -81,7 +86,11 @@ def generate_plot_grid( 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} + return {key: np.array(nr_charge[key]) for key in nr_charge} + +def compile_nr_light_data( + data_dir: str="data/", +): nr_light = { "Dataset": [], "energy": [], @@ -94,7 +103,7 @@ def generate_plot_grid( "yield_sl": [], "yield_sh": [] } - with open("data/nr_light_alt.csv","r") as file: + with open(f"{data_dir}nr_light_alt.csv","r") as file: reader = csv.reader(file, delimiter=",") next(reader) for row in reader: @@ -108,7 +117,11 @@ def generate_plot_grid( 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} + return {key: np.array(nr_light[key]) for key in nr_light} + +def compile_er_charge_data( + data_dir: str="data/", +): er_charge = { "Dataset": [], "energy": [], @@ -121,7 +134,7 @@ def generate_plot_grid( "yield_sl": [], "yield_sh": [] } - with open("data/er_charge.csv","r") as file: + with open(f"{data_dir}er_charge.csv","r") as file: reader = csv.reader(file, delimiter=",") next(reader) for row in reader: @@ -135,7 +148,11 @@ def generate_plot_grid( 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} + return {key: np.array(er_charge[key]) for key in er_charge} + +def compile_er_light_data( + data_dir: str="data/", +): er_light = { "Dataset": [], "energy": [], @@ -148,7 +165,7 @@ def generate_plot_grid( "yield_sl": [], "yield_sh": [] } - with open("data/er_light.csv","r") as file: + with open(f"{data_dir}er_light.csv","r") as file: reader = csv.reader(file, delimiter=",") next(reader) for row in reader: @@ -162,13 +179,20 @@ def generate_plot_grid( 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} + return {key: np.array(er_light[key]) for key in er_light} +def compile_lar_data( + data_dir: str="data/", +): 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 + f"{data_dir}lar_data.npz", + nr_total=compile_nr_total_data(data_dir), + nr_charge=compile_nr_charge_data(data_dir), + nr_light=compile_nr_light_data(data_dir), + er_charge=compile_er_charge_data(data_dir), + er_light=compile_er_light_data(data_dir), + ) + +if __name__ == "__main__": + + compile_lar_data() \ No newline at end of file