Merge branch 'cosmics' into 'master' (cde27218) · Commits · XRFLab109 / Allpix Squared

doc/usermanual/references.bib

+25 −0

Original line number	Diff line number	Diff line
		@@ -571,3 +571,28 @@ year = {2008}
		month = 10,
		url = {https://inspirehep.net/literature/1813192}
		}

		@INPROCEEDINGS{cry,
		author={Hagmann, Chris and Lange, David and Wright, Douglas},
		booktitle={2007 IEEE Nuclear Science Symposium Conference Record},
		title={Cosmic-ray shower generator (CRY) for Monte Carlo transport codes},
		year={2007},
		volume={2},
		number={},
		pages={1143-1146},
		doi={10.1109/NSSMIC.2007.4437209}
		}

		@online{cryphysics,
		author = {Chris Hagmann and David Lange and Doug Wright},
		title = {Monte Carlo Simulation of Proton-induced Cosmic-ray Cascades in the Atmosphere},
		subtitle = {UCRL-TM-229452},
		url = {https://nuclear.llnl.gov/simulation/doc_cry_v1.7/cry_physics.pdf}
		}

		@online{crymanual,
		author = {Chris Hagmann and David Lange and Jerome Verbeke and Doug Wright},
		title = {Cosmic-ray Shower Library (CRY)},
		subtitle = {UCRL-TM-229453},
		url = {https://nuclear.llnl.gov/simulation/doc_cry_v1.7/cry.pdf}
		}

src/modules/DepositionCosmics/CMakeLists.txt

0 → 100644

+110 −0

Original line number	Diff line number	Diff line
		# Define module and return the generated name as MODULE_NAME
		ALLPIX_UNIQUE_MODULE(MODULE_NAME)

		# The tabulated data is downloaded at compile time since the files are too big for the repository
		INCLUDE(ExternalData)
		# cmake-lint: disable=C0103
		# Define resources for external data, copy from EOS directly if available
		SET(ExternalData_URL_TEMPLATES "file:///eos/project-a/allpix-squared/www/data/%(algo)/%(hash)"
		"https://project-allpix-squared.web.cern.ch/data/%(algo)/%(hash)")

		# Make copies, not symlinks so that the files can be picked up and installed by the install target
		SET(ExternalData_NO_SYMLINKS 1)

		# Register the data files via placeholders with their name, the extension .md5 and their MD5 has as content:
		EXTERNALDATA_EXPAND_ARGUMENTS(fetch_cry_data _data_file DATA{data/cosmics_0.data})
		EXTERNALDATA_EXPAND_ARGUMENTS(fetch_cry_data _data_file DATA{data/cosmics_2100.data})
		EXTERNALDATA_EXPAND_ARGUMENTS(fetch_cry_data _data_file DATA{data/cosmics_11300.data})

		# Set up a CMake target for downloading the files, exclude from "make all" but make this module target depend on it:
		EXTERNALDATA_ADD_TARGET(fetch_cry_data)
		SET_TARGET_PROPERTIES(fetch_cry_data PROPERTIES EXCLUDE_FROM_ALL TRUE)
		ADD_DEPENDENCIES(${MODULE_NAME} fetch_cry_data)

		# Prefer local data files if not installed to other location
		GET_FILENAME_COMPONENT(ABSOLUTE_INSTALL_PREFIX ${CMAKE_INSTALL_PREFIX} ABSOLUTE)
		IF(ABSOLUTE_INSTALL_PREFIX STREQUAL PROJECT_SOURCE_DIR)
		# Use local data directory (unless other is given)
		SET(CRY_DATA_DIRECTORY
		"${CMAKE_CURRENT_BINARY_DIR}/data"
		CACHE PATH "CRY data directory" FORCE)
		ELSE()
		# Set the install path to share
		SET(CRY_DATA_DIRECTORY
		"share/${CMAKE_PROJECT_NAME}/data"
		CACHE PATH "CRY data directory" FORCE)
		ENDIF()

		SET(DATA_INSTALL_DIRECTORY "${CRY_DATA_DIRECTORY}")
		IF(NOT IS_ABSOLUTE ${CRY_DATA_DIRECTORY})
		SET(CRY_DATA_DIRECTORY
		"${CMAKE_INSTALL_PREFIX}/${CRY_DATA_DIRECTORY}"
		CACHE PATH "CRY data directory" FORCE)
		ENDIF()

		# Install the files to the data directory
		# NOTE: With default install path this does not change anything
		INSTALL(DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/data/ DESTINATION ${CRY_DATA_DIRECTORY}
		FILES_MATCHING PATTERN "*.data")

		ADD_DEFINITIONS(-DALLPIX_CRY_DATA_DIRECTORY="${CRY_DATA_DIRECTORY}")


		# Required for CMake when Geant4 has HDF5 enabled
		ENABLE_LANGUAGE(C)

		# Geant4 is required for geometry description and charge deposition.
		FIND_PACKAGE(Geant4)
		IF(NOT Geant4_FOUND)
		MESSAGE(FATAL_ERROR "Could not find Geant4, make sure to source the Geant4 environment\n"
		"$ source YOUR_GEANT4_DIR/bin/geant4.sh")
		ENDIF()

		# Add "geant4.sh" as runtime dependency for setup.sh file:
		ADD_RUNTIME_DEP(geant4.sh)

		# Add Geant4 flags before our own flags
		ADD_DEFINITIONS(${Geant4_DEFINITIONS})
		SET(CMAKE_CXX_FLAGS "${Geant4_CXX_FLAGS} ${CMAKE_CXX_FLAGS}")
		IF(CMAKE_BUILD_TYPE MATCHES DEBUG)
		SET(CMAKE_CXX_FLAGS "${Geant4_CXX_FLAGS_DEBUG} ${CMAKE_CXX_FLAGS}")
		ELSEIF(CMAKE_BUILD_TYPE MATCHES RELEASE)
		SET(CMAKE_CXX_FLAGS "${Geant4_CXX_FLAGS_RELEASE} ${CMAKE_CXX_FLAGS}")
		ENDIF()

		# Add source files to library
		ALLPIX_MODULE_SOURCES(${MODULE_NAME}
		DepositionCosmicsModule.cpp
		CosmicsGeneratorActionG4.cpp
		RNGWrapper.cpp
		cry/CRYAbsFunction.cc
		cry/CRYAbsParameter.cc
		cry/CRYBinning.cc
		cry/CRYCosLatitudeFunction.cc
		cry/CRYData.cc
		cry/CRYFunctionDict.cc
		cry/CRYGenerator.cc
		cry/CRYParameter.cc
		cry/CRYParamI.cc
		cry/CRYParticle.cc
		cry/CRYPdf.cc
		cry/CRYPrimary.cc
		cry/CRYPrimarySpectrumFunction.cc
		cry/CRYSetup.cc
		cry/CRYUtils.cc
		cry/CRYWeightFunc.cc)

		# Include Geant4 directories (NOTE Geant4_USE_FILE is not used!)
		TARGET_INCLUDE_DIRECTORIES(${MODULE_NAME} SYSTEM PRIVATE ${Geant4_INCLUDE_DIRS})

		# Add Geant4 libraries
		TARGET_LINK_LIBRARIES(${MODULE_NAME} ${Geant4_LIBRARIES} AllpixModuleDepositionGeant4)

		# Allpix Geant4 interface is required for this module
		ALLPIX_MODULE_REQUIRE_GEANT4_INTERFACE(${MODULE_NAME})

		# Register module tests
		ALLPIX_MODULE_TESTS(${MODULE_NAME} "tests")

		# Provide standard install target
		ALLPIX_MODULE_INSTALL(${MODULE_NAME})

src/modules/DepositionCosmics/CosmicsGeneratorActionG4.cpp

0 → 100644

+88 −0

Original line number	Diff line number	Diff line
		/**
		* @file
		* @brief Implements the interface of the Geant4 ParticleGun to CRY
		* @copyright Copyright (c) 2021 CERN and the Allpix Squared authors.
		* This software is distributed under the terms of the MIT License, copied verbatim in the file "LICENSE.md".
		* In applying this license, CERN does not waive the privileges and immunities granted to it by virtue of its status as an
		* Intergovernmental Organization or submit itself to any jurisdiction.
		*/

		#include "CosmicsGeneratorActionG4.hpp"
		#include "DepositionCosmicsModule.hpp"
		#include "RNGWrapper.hpp"

		#include <limits>
		#include <memory>
		#include <regex>

		#include <G4Event.hh>
		#include <G4ParticleTable.hh>

		#include "core/config/exceptions.h"
		#include "core/utils/log.h"
		#include "tools/geant4/geant4.h"

		using namespace allpix;

		CosmicsGeneratorActionG4::CosmicsGeneratorActionG4(const Configuration& config)
		: particle_gun_(std::make_unique<G4ParticleGun>()), config_(config) {

		LOG(DEBUG) << "Setting up CRY generator";
		LOG(DEBUG) << "CRY configuration: " << config_.get<std::string>("_cry_config");
		LOG(DEBUG) << "CRY data: " << config_.get<std::string>("data_path");
		auto* setup = new CRYSetup(config_.get<std::string>("_cry_config"), config_.get<std::string>("data_path"));
		cry_generator_ = std::make_unique<CRYGenerator>(setup);

		// Set up random number generator to use the Geant4-internal one which is seeded per-event:
		LOG(DEBUG) << "Configuring CRY random engine to use Geant4's event-seeded engine";
		RNGWrapper<CLHEP::HepRandomEngine>::set(CLHEP::HepRandom::getTheEngine(), &CLHEP::HepRandomEngine::flat);
		setup->setRandomFunction(RNGWrapper<CLHEP::HepRandomEngine>::rng);

		// Parse other configuration parameters:
		reset_particle_time_ = config_.get<bool>("reset_particle_time");
		}

		/**
		* Called automatically for every event
		*/
		void CosmicsGeneratorActionG4::GeneratePrimaries(G4Event* event) {

		// Let CRY generate the particles
		std::vector<CRYParticle*> vect;
		LOG(DEBUG) << "Absolute time simulated before shower: "
		<< Units::display(Units::get(cry_generator_->timeSimulated(), "s"), {"ns", "us", "ms"});
		cry_generator_->genEvent(&vect);
		LOG(DEBUG) << "CRY generated " << vect.size() << " particles";
		LOG(INFO) << "Absolute time simulated by CRY after shower: "
		<< Units::display(Units::get(cry_generator_->timeSimulated(), "s"), {"ns", "us", "ms"});

		// Update simulation time in the framework base units
		DepositionCosmicsModule::cry_instance_time_simulated_ = cry_generator_->timeSimulated() * 1e9;

		// Event time frame starts with first particle arriving
		double event_starting_time = std::numeric_limits<double>::max();
		if(!reset_particle_time_) {
		for(const auto& particle : vect) {
		event_starting_time = std::min(event_starting_time, particle->t());
		}
		}

		for(const auto& particle : vect) {
		auto* pdg_table = G4ParticleTable::GetParticleTable();
		particle_gun_->SetParticleDefinition(pdg_table->FindParticle(particle->PDGid()));
		particle_gun_->SetParticleEnergy(particle->ke() * CLHEP::MeV);
		particle_gun_->SetParticlePosition(
		G4ThreeVector(particle->x() * CLHEP::m, particle->y() * CLHEP::m, particle->z() * CLHEP::m));
		particle_gun_->SetParticleMomentumDirection(G4ThreeVector(particle->u(), particle->v(), particle->w()));

		double time = (reset_particle_time_ ? 0. : particle->t() - event_starting_time);
		particle_gun_->SetParticleTime(time);
		particle_gun_->GeneratePrimaryVertex(event);

		LOG(DEBUG) << " " << CRYUtils::partName(particle->id()) << ": charge=" << particle->charge() << std::setprecision(4)
		<< " energy=" << Units::display(particle->ke(), {"MeV", "GeV"}) << " pos="
		<< Units::display(G4ThreeVector(particle->x() * 1e3, particle->y() * 1e3, particle->z() * 1e3), {"m"})
		<< " dir. cos=" << G4ThreeVector(particle->u(), particle->v(), particle->w())
		<< " t=" << Units::display(Units::get(time, "s"), {"ns", "us", "ms"});
		}
		}

src/modules/DepositionCosmics/CosmicsGeneratorActionG4.hpp

0 → 100644

+58 −0

Original line number	Diff line number	Diff line
		/**
		* @file
		* @brief Defines the CRY interface to Geant4
		* @copyright Copyright (c) 2021 CERN and the Allpix Squared authors.
		* This software is distributed under the terms of the MIT License, copied verbatim in the file "LICENSE.md".
		* In applying this license, CERN does not waive the privileges and immunities granted to it by virtue of its status as an
		* Intergovernmental Organization or submit itself to any jurisdiction.
		*/

		#ifndef ALLPIX_COSMICS_DEPOSITION_MODULE_GENERATOR_ACTION_H
		#define ALLPIX_COSMICS_DEPOSITION_MODULE_GENERATOR_ACTION_H

		#include <memory>

		#include <G4DataVector.hh>
		#include <G4ParticleGun.hh>
		#include <G4ThreeVector.hh>
		#include <G4VUserPrimaryGeneratorAction.hh>

		#include "cry/CRYGenerator.h"
		#include "cry/CRYParticle.h"
		#include "cry/CRYSetup.h"
		#include "cry/CRYUtils.h"

		#include "core/config/Configuration.hpp"

		namespace allpix {
		/**
		* @brief Generates the particles in every event
		*/
		class CosmicsGeneratorActionG4 : public G4VUserPrimaryGeneratorAction {
		public:
		/**
		* @brief Constructs the generator action
		* @param config Configuration of the \ref DepositionCosmicsModule module
		*/
		explicit CosmicsGeneratorActionG4(const Configuration& config);

		/**
		* @brief Generate the particle for every event
		*/
		void GeneratePrimaries(G4Event*) override;

		private:
		std::unique_ptr<G4ParticleGun> particle_gun_;
		std::unique_ptr<CRYGenerator> cry_generator_;

		bool reset_particle_time_{};
		const Configuration& config_;
		};

		class GeneratorActionInitializationMaster {
		public:
		explicit GeneratorActionInitializationMaster(const Configuration&){};
		};
		} // namespace allpix

		#endif /* ALLPIX_COSMICS_DEPOSITION_MODULE_GENERATOR_ACTION_H */

src/modules/DepositionCosmics/DepositionCosmicsModule.cpp

0 → 100644

+177 −0

Original line number	Diff line number	Diff line
		/**
		* @file
		* @brief Implementation of DepositionCosmics module
		* @copyright Copyright (c) 2021 CERN and the Allpix Squared authors.
		* This software is distributed under the terms of the MIT License, copied verbatim in the file "LICENSE.md".
		* In applying this license, CERN does not waive the privileges and immunities granted to it by virtue of its status as an
		* Intergovernmental Organization or submit itself to any jurisdiction.
		*/

		#include "DepositionCosmicsModule.hpp"
		#include "CosmicsGeneratorActionG4.hpp"

		#include "tools/geant4/MTRunManager.hpp"
		#include "tools/geant4/RunManager.hpp"
		#include "tools/geant4/geant4.h"

		#include "../DepositionGeant4/ActionInitializationG4.hpp"

		#include <filesystem>
		#include <sstream>
		#include <string>
		#include <utility>

		#include "core/utils/log.h"

		using namespace allpix;

		thread_local double DepositionCosmicsModule::cry_instance_time_simulated_ = 0;

		DepositionCosmicsModule::DepositionCosmicsModule(Configuration& config, Messenger* messenger, GeometryManager* geo_manager)
		: DepositionGeant4Module(config, messenger, geo_manager) {

		// Enable multithreading of this module if multithreading is enabled
		allow_multithreading();

		// Set defaults for configuration parameters:
		config_.setDefault("return_neutrons", true);
		config_.setDefault("return_protons", true);
		config_.setDefault("return_gammas", true);
		config_.setDefault("return_electrons", true);
		config_.setDefault("return_muons", true);
		config_.setDefault("return_pions", true);
		config_.setDefault("return_kaons", true);
		config_.setDefault("altitude", Units::get(0, "m"));
		config_.setDefault("min_particles", 1);
		config_.setDefault("max_particles", 1000000);
		config_.setDefault("latitude", 53.0);
		config_.setDefault("date", "12-31-2020");
		config_.setDefault("reset_particle_time", false);

		// Force source type and position:
		config_.set("source_type", "cosmics");
		config_.set("source_position", ROOT::Math::XYZPoint());

		// Add the particle source position to the geometry
		geo_manager_->addPoint(config_.get<ROOT::Math::XYZPoint>("source_position", ROOT::Math::XYZPoint()));

		// Register lookup path for CRY data files:
		if(config_.has("data_path")) {
		auto path = config_.getPath("data_path", true);
		if(!std::filesystem::is_directory(path)) {
		throw InvalidValueError(config_, "data_path", "path does not point to a directory");
		}
		LOG(TRACE) << "Registered CRY data path from configuration: " << path;
		} else {
		if(std::filesystem::is_directory(ALLPIX_CRY_DATA_DIRECTORY)) {
		config_.set<std::string>("data_path", ALLPIX_CRY_DATA_DIRECTORY);
		LOG(TRACE) << "Registered CRY data path from system: " << ALLPIX_CRY_DATA_DIRECTORY;
		} else {
		const char* data_dirs_env = std::getenv("XDG_DATA_DIRS");
		if(data_dirs_env == nullptr \|\| strlen(data_dirs_env) == 0) {
		data_dirs_env = "/usr/local/share/:/usr/share/:";
		}

		std::vector<std::string> data_dirs = split<std::string>(data_dirs_env, ":");
		for(auto data_dir : data_dirs) {
		if(data_dir.back() != '/') {
		data_dir += "/";
		}
		data_dir += std::string(ALLPIX_PROJECT_NAME) + std::string("/data");

		if(std::filesystem::is_directory(data_dir)) {
		config_.set<std::string>("data_path", data_dir);
		LOG(TRACE) << "Registered CRY data path from XDG_DATA_DIRS: " << data_dir;
		} else {
		throw ModuleError("Cannot find CRY data files, provide them in the configuration, via XDG_DATA_DIRS or "
		"in system directory " +
		std::string(ALLPIX_CRY_DATA_DIRECTORY));
		}
		}
		}
		}

		// Build configuration for CRY - it expects a single string with all tokens separated by whitespace:
		std::stringstream cry_config;

		cry_config << " returnNeutrons " << static_cast<int>(config_.get<bool>("return_neutrons")) << " returnProtons "
		<< static_cast<int>(config_.get<bool>("return_protons")) << " returnGammas "
		<< static_cast<int>(config_.get<bool>("return_gammas")) << " returnElectrons "
		<< static_cast<int>(config_.get<bool>("return_electrons")) << " returnMuons "
		<< static_cast<int>(config_.get<bool>("return_muons")) << " returnPions "
		<< static_cast<int>(config_.get<bool>("return_pions")) << " returnKaons "
		<< static_cast<int>(config_.get<bool>("return_kaons"));

		// Select altitude:
		auto altitude = config_.get<int>("altitude");
		if(altitude != 0 && altitude != 2100000 && altitude != 11300000) {
		throw InvalidValueError(config_, "altitude", "only altitudes of 0m, 2100m and 11300m are supported");
		}
		cry_config << " altitude " << static_cast<int>(Units::convert(altitude, "m"));

		cry_config << " nParticlesMin " << config_.get<unsigned int>("min_particles") << " nParticlesMax "
		<< config_.get<unsigned int>("max_particles");

		if(!config_.has("area")) {
		// Calculate subbox length required from the maximum coordinates of the setup:
		LOG(DEBUG) << "Calculating subbox length from setup size";
		auto min = geo_manager_->getMinimumCoordinate();
		auto max = geo_manager_->getMaximumCoordinate();
		auto size = std::max(max.x() - min.x(), max.y() - min.y());
		auto size_meters = static_cast<unsigned int>(std::ceil(Units::convert(size, "m")));
		if(size_meters > 300) {
		throw ModuleError("Size of the setup too large, tabulated data only available for areas up to 300m");
		}
		LOG(DEBUG) << "Maximum side length (in x,y): " << Units::display(size, {"mm", "cm", "m"})
		<< ", selecting subbox of size " << size_meters << "m";
		cry_config << " subboxLength " << size_meters;
		} else {
		auto area = Units::convert(config_.get<int>("area"), "m");
		if(area > 300) {
		throw InvalidValueError(config_, "area", "only areas with side lengths of up to 300m are supported");
		}
		LOG(DEBUG) << "Configuring subbox of size " << area << "m from configuration parameter";
		cry_config << " subboxLength " << area;
		}

		auto latitude = config_.get<double>("latitude");
		if(latitude < -90 \|\| latitude > 90) {
		throw InvalidValueError(config_, "latitude", "latitude has to be between 90.0 (north pole) and -90.0 (south pole)");
		}
		cry_config << " latitude " << latitude;
		cry_config << " date " << config_.get<std::string>("date");

		// Set CRY configuration string as internal key:
		config_.set<std::string>("_cry_config", cry_config.str());
		}

		void DepositionCosmicsModule::initialize_g4_action() {
		auto* action_initialization =
		new ActionInitializationG4<CosmicsGeneratorActionG4, GeneratorActionInitializationMaster>(config_);
		run_manager_g4_->SetUserInitialization(action_initialization);
		}

		void DepositionCosmicsModule::finalizeThread() {
		{
		LOG(DEBUG) << "CRY instance reports simulation time of "
		<< Units::display(cry_instance_time_simulated_, {"us", "ms", "s"});
		std::lock_guard<std::mutex> lock{stats_mutex_};
		total_time_simulated_ += cry_instance_time_simulated_;
		}

		// Call base class thread finalization:
		DepositionGeant4Module::finalizeThread();
		}

		void DepositionCosmicsModule::finalize() {

		// Without multithreading we need to fetch the total simulation time from the main thread:
		if(!multithreadingEnabled()) {
		total_time_simulated_ = cry_instance_time_simulated_;
		}

		LOG(STATUS) << "Total simulated time in CRY: " << Units::display(total_time_simulated_, {"us", "ms", "s"});

		// Call base class finalization:
		DepositionGeant4Module::finalize();
		}