Skip to content

Commit

Permalink
LArNEST fluctuations. (#165)
Browse files Browse the repository at this point in the history 
* removed calls to 'clamp' and fixed compiler warnings for non-void functions returning void

* working on binomial fluctuations for lar

* Added preliminary fluctuations model for LAr, just does binomial fluctuations the same way as is done for LXe

* fixed codefactor errors

* fixed codefactor errors
  • Loading branch information
@infophysics
infophysics authored Sep 8, 2022
1 parent fec16c1 commit 489cb24
Show file tree
Hide file tree
Showing 10 changed files with 526 additions and 79 deletions.
10 changes: 8 additions & 2 deletions examples/LArNEST/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
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)
173 changes: 173 additions & 0 deletions examples/LArNEST/LArNESTFluctuationBenchmarks.cpp
View file
Original file line number Diff line number Diff line change
@@ -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 <fstream>
#include <vector>
#include <string>
#include <iostream>
#include <stdlib.h>

#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<double> 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<NEST::LArInteraction> particle_types = {
NEST::LArInteraction::NR, NEST::LArInteraction::ER, NEST::LArInteraction::Alpha
};
std::vector<std::string> 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<double> 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<double> Nph(num_events);
std::vector<double> Ne(num_events);
std::vector<double> Nex(num_events);
std::vector<double> 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;
}
6 changes: 3 additions & 3 deletions examples/LArNEST/LArNESTBenchmarks.cpp → ...s/LArNEST/LArNESTMeanYieldsBenchmarks.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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++)
Expand Down Expand Up @@ -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";
}
}
}
Expand Down
84 changes: 84 additions & 0 deletions examples/LArNEST/LArNESTNeutronCapture.cpp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,84 @@
/**
* @file LArNESTNeutronCapture.cpp
* @author Nicholas Carrara [nmcarrara@ucdavis.edu]
* @brief
* @version 0.1
* @date 2022-05-23
*/
#include <fstream>
#include <vector>
#include <string>
#include <iostream>
#include <stdlib.h>

#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<double> electric_field = {
1, 100, 200, 500, 600, 1000, 1500, 2500, 6000, 9000, 9500
};

std::vector<double> 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;
}
2 changes: 1 addition & 1 deletion examples/LArNEST/LegacyLArNESTBenchmarks.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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++)
Expand Down
44 changes: 30 additions & 14 deletions include/NEST/LArNEST.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,11 @@ namespace NEST
Alpha = 2,
};

enum class LArFluctuationModel
{
Default = 0,
};

struct LArNRYieldsParameters
{
double alpha = {11.10};
Expand Down Expand Up @@ -175,9 +180,6 @@ namespace NEST

struct LArYieldFluctuationResult
{
double TotalYieldFluctuation;
double QuantaYieldFluctuation;
double LightYieldFluctuation;
double NphFluctuation;
double NeFluctuation;
double NexFluctuation;
Expand Down Expand Up @@ -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; }
Expand Down Expand Up @@ -275,7 +277,8 @@ namespace NEST
*
*/
LArYieldFluctuationResult GetYieldFluctuations(
const YieldResult &yields, double density
const LArYieldResult &yields,
double density
);
/**
* @brief
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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;
Expand All @@ -526,5 +540,7 @@ namespace NEST

ThomasImelParameters fThomasImelParameters;
DriftParameters fDriftParameters;

enum LArFluctuationModel fLArFluctuationModel;
};
}
Loading

0 comments on commit 489cb24

Please sign in to comment.