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 <[email protected]> Co-authored-by: Nicholas Carrara <[email protected]>
NESTCollaboration · May 5, 2022 · 85133c3 · 85133c3
1 parent f1009d7
commit 85133c3
Show file tree

Hide file tree

Showing 17 changed files with 2,203 additions and 2 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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)
diff --git a/examples/LArNEST/LArNESTBenchmarks.cpp b/examples/LArNEST/LArNESTBenchmarks.cpp
@@ -0,0 +1,113 @@
+/**
+ * @file legacyLArNEST.cpp
+ * @author NEST Collaboration
+ * @author Nicholas Carrara [[email protected]]
+ * @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;
+    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<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("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<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
+            });
+        }
+        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;
+}
diff --git a/examples/LArNEST/LegacyLArNESTBenchmarks.cpp b/examples/LArNEST/LegacyLArNESTBenchmarks.cpp
@@ -0,0 +1,113 @@
+/**
+ * @file legacyLArNEST.cpp
+ * @author NEST Collaboration
+ * @author Nicholas Carrara [[email protected]]
+ * @brief 
+ * @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 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<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);
+    }
+
+    int pdgcode = atoi(argv[2]);
+    double density = atof(argv[3]);
+    std::vector<double> 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<double> TotalYield(num_events);
+            std::vector<double> QuantaYield(num_events);
+            std::vector<double> LightYield(num_events);
+            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.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;
+}
diff --git a/include/Detectors/LArDetector.hh b/include/Detectors/LArDetector.hh
@@ -0,0 +1,24 @@
+/**
+ * @file LArDetector.hh
+ * @author NEST Collaboration
+ * @author Nicholas Carrara [[email protected]]
+ * @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:
+};