Lets you read in and simulate 2D shower profiles generated from Corsika

examples/CMakeLists.txt

@@ -1,5 +1,6 @@
 cmake_minimum_required( VERSION 3.20 FATAL_ERROR )
 add_subdirectory(shower_profile)
+add_subdirectory(shower_profile_2d)
 add_subdirectory(electric_dipole_antenna)
 add_subdirectory(coulomb)
 add_subdirectory(point_charge_coulomb)

examples/shower_profile/01generate_weighting_field.cpp

@@ -19,15 +19,15 @@ int main(int argc, char* argv[]) {
   }
 
   // Domain of weighting field
-  CoordVector start_coords = CU::MakeCoordVectorTRZ(-500.0, -10.0, -551.0);
-  CoordVector end_coords = CU::MakeCoordVectorTRZ(1350.0, 700.0, -549.0);
+  CoordVector start_coords = CU::MakeCoordVectorTRZ(-500.0, -10.0, -554.);
+  CoordVector end_coords = CU::MakeCoordVectorTRZ(1350.0, 700.0, -546.);
 
   // Filter parameters
   scalar_t tp = 1.0;
   unsigned int N = 6;
 
   // Sampling parameters
-  scalar_t os_factor = 5;
+  scalar_t os_factor = 3;
   scalar_t r_min = 0.1;
 
   scalar_t index_of_refraction = 1.78;

examples/shower_profile_2d/02calculate_shower_emission2d.cpp

@@ -0,0 +1,66 @@
+#include <iostream>
+#include <fstream>
+#include <algorithm>
+#include <vector>
+#include <stdlib.h>
+#include "Eisvogel/Common.hh"
+#include "Eisvogel/SignalCalculator.hh"
+#include "Eisvogel/Current0D.hh"
+#include "Eisvogel/SignalExport.hh"
+#include "Eisvogel/CoordUtils.hh"
+#include "Eisvogel/SparseCurrentDensity3D.hh"
+#include "shower_creator.h"
+#include "shower_2D.h"
+#include "constants.h"
+#include "units.h"
+
+using CurrentElement = std::pair<CoordVector, FieldVector>;
+
+
+int main (int argc, char* argv[]) {
+    scalar_t shower_energy = 1.e18;                  // energy of the shower, in eV
+    scalar_t shower_zenith = 90 * units::degree;     // zenith angle of the shower. I recommend leaving this as is
+    scalar_t shower_azimuth = 0;                     // azimuth of the shower. I also recommend leaving this
+    int is_hadronic = 0;                             // sets the type of the shower. Choose 1 for hardonic or 0 for electromagnetic shower
+    int i_shower = 9;                                // index of the shower profile to be chosen from the library. Choose a value between 0 and 9 or set to NULL to have a shower selected at random
+    std::array<float, 3> shower_vertex = {-106  , .1, -165};    // position of the shower vertex. If you change this, make sure to also adjust the weighting field accordingly
+    scalar_t sampling_rate = 5.;                     // sampling rate that is used when calculating the radio signal
+
+    std::string wf_path;
+    std::string output_path;
+    if (argc < 2) {
+        wf_path = "weighting_field";
+    } else {
+        wf_path = argv[1];
+    }
+    if (argc < 3) {
+        output_path = "shower_radio_emission.csv";
+    } else {
+        output_path = argv[2];
+    }
+    SignalCalculator signal_calc(wf_path);
+    std::cout << "Done reading weighting field \n";
+
+    showers::ShowerCreator2D shower_creator("/home/welling/RadioNeutrino/scripts/corsika/plot_examples/charge_excess_profile.hdf5");
+    showers::Shower2D shower = shower_creator.create_shower(shower_vertex, shower_energy, shower_zenith, shower_azimuth, is_hadronic);
+
+    std::vector<scalar_t> signal_times, signal_values;
+    DeltaVector delta_vec =  CoordUtils::MakeDeltaVectorTXYZ(.5, .5, .5, .5);
+/*
+    std::cout << "Getting current \n";
+    SparseCurrentDensity3D current = shower.get_current(delta_vec);
+    std::cout << "Starting radio signal calculation \n";
+    for(scalar_t cur_t = 1050; cur_t < 1300; cur_t += 1. / sampling_rate) {
+        // std::cout << cur_t << "\n";
+        scalar_t cur_signal = signal_calc.ComputeSignal(current, cur_t);
+        signal_times.push_back(cur_t);
+        signal_values.push_back(cur_signal);
+    }
+
+    for (int i; i < signal_values.size(); i++) {
+        signal_values[i] = signal_values[i] * (1. / constants::epsilon_0 / constants::c / constants::c);
+    }
+    ExportSignal(signal_times, signal_values, output_path);
+    */
+    return 0;
+}

examples/shower_profile_2d/CMakeLists.txt

@@ -0,0 +1,13 @@
+set(CMAKE_CXX_STANDARD 20)
+set(CMAKE_CXX_STANDARD_REQUIRED YES)
+
+add_executable(02shower_simulation2d 02calculate_shower_emission2d.cpp)
+target_link_libraries(02shower_simulation2d eisvogel)
+target_link_libraries(02shower_simulation2d ShowerCreator)
+
+
+target_include_directories(02shower_simulation2d PUBLIC
+                          "${PROJECT_SOURCE_DIR}/extern/shower_profile/shower_creator"
+                          )
+configure_file(${PROJECT_SOURCE_DIR}/examples/shower_profile/plot_results.py ${CMAKE_CURRENT_BINARY_DIR} COPYONLY
+)

extern/shower_profile/CMakeLists.txt

@@ -7,12 +7,10 @@ add_library(Shower2D SHARED shower_2D/shower_2D.cpp)
 add_library(Environment SHARED environment/ice_profile.cpp)
 add_library(ShowerCreator SHARED shower_creator/shower_creator.cpp)
 
-target_link_libraries(ShowerCreator Shower1D)
-target_link_libraries(ShowerCreator Shower2D)
-target_link_libraries(ShowerCreator Environment)
+find_package(HDF5 REQUIRED)
+target_link_libraries(ShowerCreator PUBLIC Shower1D Shower2D Environment ${HDF5_LIBRARY_DIRS} hdf5_serial_cpp)
 target_link_libraries(Shower1D PUBLIC Environment)
 target_link_libraries(Shower2D PUBLIC Environment)
-
 target_include_directories(Shower1D PUBLIC
                           "${PROJECT_SOURCE_DIR}/extern/shower_profile/environment"
                           "${PROJECT_SOURCE_DIR}/extern/shower_profile/charge_excess_profile"
@@ -25,6 +23,7 @@ target_include_directories(ShowerCreator PUBLIC
                           "${PROJECT_SOURCE_DIR}/extern/shower_profile/shower_1D"
                           "${PROJECT_SOURCE_DIR}/extern/shower_profile/shower_2D"
                           "${PROJECT_SOURCE_DIR}/extern/utilities"
+                          "${HDF5_INCLUDE_DIRS}"
                           )
 target_include_directories(Environment PUBLIC
                           "${PROJECT_SOURCE_DIR}/extern/utilities"

extern/shower_profile/charge_excess_profile/charge_excess_profile.cpp

@@ -27,16 +27,17 @@ class ChargeExcessProfile2D {
 
 	double get_charge_excess(double gram, double rad) {
 		int radius_index = 0;
+		//std::cout << "---> " << gram << ", " << grammage[2] << "\n";
 		for (int i_radius=1;i_radius < radius.size(); i_radius++) {
 			radius_index = i_radius - 1;
-			if (rad > radius[i_radius]) {
+			if (rad < radius[i_radius]) {
 				break;
 			}
 		}
 		int grammage_index = 0;
 		for (int i_grammage = 1; i_grammage < grammage.size() - 1; i_grammage++) {
 			grammage_index = i_grammage - 1;
-			if (gram > grammage[i_grammage]) {
+			if (gram < grammage[i_grammage]) {
 				break;
 			}
 		}
@@ -46,8 +47,8 @@ class ChargeExcessProfile2D {
 
 		double delta_ce_2 = charge_excess(grammage_index + 1, radius_index + 1) - charge_excess(grammage_index, radius_index + 1);
 		double grammage_interpolation_2 = charge_excess(grammage_index, radius_index + 1) + delta_ce_2 / delta_grammage * (gram - grammage[grammage_index]);
-
 		double delta_r = radius[radius_index + 1] - radius[radius_index];
+
 		return grammage_interpolation_1 + (grammage_interpolation_2 - grammage_interpolation_1) / delta_r * (rad - radius[radius_index]);
 	}
 };

extern/shower_profile/convert_shower_profile.py

@@ -51,7 +51,7 @@ def write_shower(f, E, type, grammage, q):
 		group_2 = group.create_group(str(key_2))
 		group_2.create_dataset('charge_excess', data=profiles)
 		group_2.create_dataset('depth', data=list(shower_library[key_1][key_2]['depth']))
-
+		print(np.max(shower_library[key_1][key_2]['depth']))
 	group.create_dataset('energies', data=list(np.log10(energies)))
 hf.close()
 json.dump(output_dict, open('ARZ_library_v1.2.json', 'w'))

extern/shower_profile/shower_2D/shower_2D.cpp

@@ -33,6 +33,7 @@ SparseCurrentDensity3D showers::Shower2D::get_current(DeltaVector voxel_size) {
     scalar_t delta_t = getT(voxel_size);
     double max_grammage = charge_excess_profile.grammage[charge_excess_profile.grammage.size() - 1];
     double max_radius = charge_excess_profile.radius[charge_excess_profile.radius.size() - 1];
+    std::cout << "Max Radius: " << max_radius << "\n";
     double integrated_grammage = 0;
     double delta_s = delta_t * constants::c;
     double shower_t = 0;
@@ -45,10 +46,9 @@ SparseCurrentDensity3D showers::Shower2D::get_current(DeltaVector voxel_size) {
         cos(zenith) * sin(azimuth),
         sin(zenith) * cos(azimuth) * cos(azimuth) - sin(zenith) * sin(azimuth) * sin(azimuth)
         };
-    std::cout << "vertical: " << axis_normal_ver[0] << ", " << axis_normal_ver[1] << ", " << axis_normal_ver[2] << "\n";
     double current_radius;
     double current_ce;
-    while (integrated_grammage < max_grammage) {
+    while (integrated_grammage < max_grammage / 2) {
         integrated_grammage = integrated_grammage + ice_profile.get_density(
             shower_pos[0], shower_pos[1], shower_pos[2]
             ) * delta_s;
@@ -59,6 +59,7 @@ SparseCurrentDensity3D showers::Shower2D::get_current(DeltaVector voxel_size) {
         for (float hor_offset = -max_radius; hor_offset < max_radius; hor_offset = hor_offset + delta_s) {
             for (float ver_offset = -max_radius; ver_offset < max_radius; ver_offset = ver_offset + delta_s) {
                 current_radius = sqrt(hor_offset * hor_offset + ver_offset * ver_offset);
+
                 if (current_radius < max_radius) {
                     current_ce = charge_excess_profile.get_charge_excess(integrated_grammage, current_radius);
                     if (current_ce > 0) {

extern/shower_profile/shower_creator/shower_creator.cpp

@@ -1,4 +1,5 @@
 #include "Eisvogel/Common.hh"
+#include "Eisvogel/NDArray.hh"
 #include "shower_creator.h"
 #include "shower_1D.h"
 #include "shower_2D.h"
@@ -10,6 +11,7 @@
 #include <time.h>
 #include <stdexcept>
 #include <random>
+#include "H5Cpp.h"
 
 
 showers::ShowerCreator::ShowerCreator(
@@ -105,67 +107,61 @@ showers::ShowerCreator2D::ShowerCreator2D(
 ) {
 	shower_file = file_path;
 	density_profile = environment::IceProfile();
-	int reading_complete = 0;
-	int N;
-	FILE * pFile = fopen(shower_file.c_str(), "rb");
-	ce_profiles[0] = std::map<double, std::vector<showers::ChargeExcessProfile2D>>();
-	ce_profiles[1] = std::map<double, std::vector<showers::ChargeExcessProfile2D>>();
-
-	while (reading_complete == 0) {
-		if(fread(&N, sizeof(uint32_t), 1, pFile)==0) {
-			reading_complete = 1;
-		} else {
-
-		showers::ChargeExcessProfile2D profile = read_shower(
-				pFile,
-				&N
-				);
-
-		int find_shower = ce_profiles[profile.hadronic].count(profile.energy);
-		if (find_shower == 0) {
-			ce_profiles[profile.hadronic][profile.energy] = std::vector<showers::ChargeExcessProfile2D>();
-			stored_energies[profile.hadronic].push_back(profile.energy);
-		}
-		ce_profiles[profile.hadronic][profile.energy].push_back(profile);
-		}
-	}
-
-	fclose(pFile);
+	ce_profiles.push_back(
+		read_shower(file_path)
+	);
+
 }
 
 showers::ChargeExcessProfile2D showers::ShowerCreator2D::read_shower(
-	FILE *f,
-	int *N
+	std::string file_path
 	) {
-    showers::ChargeExcessProfile2D profile({*N, 10});
-    profile.radius.resize(10);  // Resize the radius vector based on the read size
-
-    for (int i = 1; i < profile.radius.size(); i++) {
-        profile.radius[i] = profile.radius[i - 1] + 2 * units::cm;
-    }
-
-
-    fread(&profile.hadronic, sizeof(uint32_t), 1, f);
-    fread(&profile.energy, sizeof(double), 1, f);
-    profile.grammage.resize(*N);
-
-    scalar_t r_max = profile.radius[profile.radius.size() - 1];
-
-    std::vector<double> charge_io;
-    charge_io.resize(*N);  // Resize the charge_io vector based on the read size
-
-    fread(&(profile.grammage[0]), sizeof(double), *N, f);
-    fread(&(charge_io[0]), sizeof(double), *N, f);
-
-    for (int i = 0; i < profile.grammage.size(); i++) {
-        profile.grammage[i] = profile.grammage[i] * units::kilogram / units::square_meter;
-        for (int j=0; j < profile.radius.size(); j ++) {
-            profile.charge_excess(i, j) = charge_io[i] * (r_max - profile.radius[j]);
-        }
-    }
+	H5::H5File this_file(file_path, H5F_ACC_RDONLY);
+    H5::DataSet this_dataset = this_file.openDataSet("ce");
+	H5T_class_t type_class = this_dataset.getTypeClass();
+	H5::DataSpace dataspace = this_dataset.getSpace();
+	int rank = dataspace.getSimpleExtentNdims();
+	hsize_t dimensions_out[2];
+	int dim = dataspace.getSimpleExtentDims(dimensions_out, NULL);
+	hsize_t start[2] = {0, 0};
+	dataspace.selectHyperslab(H5S_SELECT_SET, dimensions_out, start);
+	H5::DataSpace memspace(rank, dimensions_out);
+	memspace.selectHyperslab(H5S_SELECT_SET, dimensions_out, start);
+	double data_out[dimensions_out[0]][dimensions_out[1]] = {0};
+	this_dataset.read( data_out, H5::PredType::NATIVE_DOUBLE, memspace, dataspace );	
+	showers::ChargeExcessProfile2D profile({dimensions_out[0], dimensions_out[1]});
+	for(int i=0; i < dimensions_out[0]; i++) {
+		for (int j=0; j < dimensions_out[1]; j++) {
+			profile.charge_excess(i, j) = data_out[i][j];
+		}
+	}
+	H5::DataSet grammage_dataset = this_file.openDataSet("depth");
+	H5::DataSet radius_dataset = this_file.openDataSet("radius");
+	profile.grammage = readDataSet(&grammage_dataset);
+	// for (int i=0;i<profile.grammage.size();i++) {
+	// 	profile.grammage[i] = profile.grammage[i] * units::kilogram / units::square_meter;
+	// }
+	profile.radius = readDataSet(&radius_dataset);
     return profile;
 }
 
+std::vector<double> showers::ShowerCreator2D::readDataSet(
+	H5::DataSet *dataset
+) {
+	H5T_class_t type_class = dataset->getTypeClass();
+	H5::DataSpace dataspace = dataset->getSpace();
+	int rank = dataspace.getSimpleExtentNdims();
+	hsize_t dimensions_out[rank];
+	int dim = dataspace.getSimpleExtentDims(dimensions_out, NULL);
+	hsize_t start[rank] = {0};
+	dataspace.selectHyperslab(H5S_SELECT_SET, dimensions_out, start);
+	H5::DataSpace memspace(rank, dimensions_out);
+	memspace.selectHyperslab(H5S_SELECT_SET, dimensions_out, start);
+	double data_out[dimensions_out[0]];
+	dataset->read(data_out, H5::PredType::NATIVE_DOUBLE, memspace, dataspace);
+	std::vector<double> vec_out (data_out, data_out + dimensions_out[0]);
+	return vec_out;
+}
 
 showers::Shower2D showers::ShowerCreator2D::create_shower(
 	std::array<float, 3> pos,
@@ -174,25 +170,25 @@ showers::Shower2D showers::ShowerCreator2D::create_shower(
 		double az,
 		int had
 ) {
-	double closest_energy = stored_energies[had][0];
-	double energy_diff = std::abs(en - stored_energies[had][0]);
-	for (int i=1; i < stored_energies[had].size(); i++) {
-		if (std::abs(en - stored_energies[had][i]) < energy_diff) {
-			closest_energy = stored_energies[had][i];
-			energy_diff = std::abs(en - stored_energies[had][i]);
-		}
-	}
-	std::default_random_engine generator{static_cast<long unsigned int>(time(NULL))};;
-	std::uniform_int_distribution<int> dist(0, ce_profiles[had][closest_energy].size());
-        int i_shower = 4;
-        std::cout << "Pick shower " << i_shower << " out of " << ce_profiles[had][closest_energy].size() << " showers. \n";
+	// double closest_energy = stored_energies[had][0];
+	// double energy_diff = std::abs(en - stored_energies[had][0]);
+	// for (int i=1; i < stored_energies[had].size(); i++) {
+	// 	if (std::abs(en - stored_energies[had][i]) < energy_diff) {
+	// 		closest_energy = stored_energies[had][i];
+	// 		energy_diff = std::abs(en - stored_energies[had][i]);
+	// 	}
+	// }
+	// std::default_random_engine generator{static_cast<long unsigned int>(time(NULL))};;
+	// std::uniform_int_distribution<int> dist(0, ce_profiles[had][closest_energy].size());
+    //     int i_shower = 4;
+    //     std::cout << "Pick shower " << i_shower << " out of " << ce_profiles[had][closest_energy].size() << " showers. \n";
 	return showers::Shower2D(
 			pos,
 			en,
 			zen,
 			az,
-			ce_profiles[had][closest_energy][i_shower],
-			en / closest_energy,
+			ce_profiles[0], //[had][closest_energy][i_shower],
+			1, //en / closest_energy,
 			density_profile
 	);
-};
+};

extern/shower_profile/shower_creator/shower_creator.h

@@ -9,6 +9,8 @@
 #include <string>
 #include <map>
 #include <random>
+#include "H5Cpp.h"
+
 
 namespace showers {
 class ShowerCreator {
@@ -48,14 +50,15 @@ class ShowerCreator2D {
 			int had
 	);
 	ChargeExcessProfile2D read_shower(
-		FILE *f,
-		int *N
+		std::string file_path
 		);
 private:
 	std::string shower_file;
 	environment::IceProfile density_profile;
-	std::map<int, std::map<double, std::vector<showers::ChargeExcessProfile2D>>>ce_profiles;
-	std::map<int, std::vector<double>> stored_energies;
+	std::vector<showers::ChargeExcessProfile2D> ce_profiles;
+	std::vector<double> readDataSet(
+		H5::DataSet *dataset
+	);
 };