Merge branch 'fix_induction_gain' into 'master' (57727089) · Commits · XRFLab109 / Allpix Squared

src/modules/GenericPropagation/GenericPropagationModule.cpp

+222 −83

Original line number	Diff line number	Diff line
		@@ -97,6 +97,7 @@ GenericPropagationModule::GenericPropagationModule(Configuration& config,
		// Set defaults for charge carrier multiplication
		config_.setDefault<std::string>("multiplication_model", "none");
		config_.setDefault<double>("multiplication_threshold", 1e-2);
		config_.setDefault<unsigned int>("max_multiplication_level", 5);

		// Copy some variables from configuration to avoid lookups:
		temperature_ = config_.get<double>("temperature");
		@@ -116,6 +117,7 @@ GenericPropagationModule::GenericPropagationModule(Configuration& config,
		propagate_holes_ = config_.get<bool>("propagate_holes");
		charge_per_step_ = config_.get<unsigned int>("charge_per_step");
		max_charge_groups_ = config_.get<unsigned int>("max_charge_groups");
		max_multiplication_level_ = config.get<unsigned int>("max_multiplication_level");

		// Enable multithreading of this module if multithreading is enabled and no per-event output plots are requested:
		// FIXME: Review if this is really the case or we can still use multithreading
		@@ -227,8 +229,44 @@ void GenericPropagationModule::initialize() {
		static_cast<int>(Units::convert(integration_time_, "ns") * 5),
		0,
		static_cast<double>(Units::convert(integration_time_, "ns")));
		gain_histo_ = CreateHistogram<TH1D>(
		"gain_histo", "Gain per charge carrier group after propagation;gain;number of groups transported", 500, 1, 25);
		if(!multiplication_.is<NoImpactIonization>()) {
		gain_primary_histo_ = CreateHistogram<TH1D>(
		"gain_primary_histo",
		"Gain per primarily induced charge carrier group after propagation;gain;number of groups transported",
		500,
		1,
		25);
		gain_all_histo_ =
		CreateHistogram<TH1D>("gain_all_histo",
		"Gain per charge carrier group after propagation;gain;number of groups transported",
		500,
		1,
		25);
		gain_e_histo_ =
		CreateHistogram<TH1D>("gain_e_histo",
		"Gain per primary electron group after propagation;gain;number of groups transported",
		500,
		1,
		25);
		gain_h_histo_ =
		CreateHistogram<TH1D>("gain_h_histo",
		"Gain per primary hole group after propagation;gain;number of groups transported",
		500,
		1,
		25);
		multiplication_level_histo_ = CreateHistogram<TH1D>(
		"multiplication_level_histo",
		"Multiplication level of propagated charge carriers;multiplication level;charge carriers",
		static_cast<int>(max_multiplication_level_),
		0,
		static_cast<int>(max_multiplication_level_));
		multiplication_depth_histo_ =
		CreateHistogram<TH1D>("multiplication_depth_histo",
		"Generation depth of charge carriers via impact ionization;depth [mm];charge carriers",
		200,
		-model_->getSensorSize().z() / 2.,
		model_->getSensorSize().z() / 2.);
		}
		}

		// Prepare mobility model
		@@ -246,6 +284,12 @@ void GenericPropagationModule::initialize() {
		<< "This might lead to unphysical gain values.";
		}

		// Check for propagating both types of charge carrier
		if(!multiplication_.is<NoImpactIonization>() && (!propagate_electrons_ \|\| !propagate_holes_)) {
		LOG(WARNING) << "Not propagating both types of charge carriers with charge multiplication enabled may lead to "
		"unphysical results";
		}

		// Prepare trapping model
		trapping_ = Trapping(config_);

		@@ -309,62 +353,24 @@ void GenericPropagationModule::run(Event* event) {
		}
		charges_remaining -= charge_per_step;

		// Get position and propagate through sensor
		auto initial_position = deposit.getLocalPosition();

		// Add point of deposition to the output plots if requested
		if(output_linegraphs_) {
		output_plot_points.emplace_back(
		std::make_tuple(deposit.getGlobalTime(), charge_per_step, deposit.getType(), CarrierState::MOTION),
		std::vector<ROOT::Math::XYZPoint>());
		}

		// Propagate a single charge deposit
		auto [final_position, time, gain, state] = propagate(initial_position,
		auto [recombined, trapped, propagated, steps, time] = propagate(event,
		deposit,
		deposit.getLocalPosition(),
		deposit.getType(),
		charge_per_step,
		deposit.getLocalTime(),
		event->getRandomEngine(),
		output_plot_points,
		charge_per_step);

		if(state == CarrierState::RECOMBINED) {
		LOG(DEBUG) << " Recombined " << charge_per_step << " at " << Units::display(final_position, {"mm", "um"})
		<< " in " << Units::display(time, "ns") << " time, removing";
		recombined_charges_count += charge_per_step;
		if(output_plots_) {
		recombination_time_histo_->Fill(static_cast<double>(Units::convert(time, "ns")), charge_per_step);
		}
		} else if(state == CarrierState::TRAPPED) {
		LOG(DEBUG) << " Trapped " << charge_per_step << " at " << Units::display(final_position, {"mm", "um"})
		<< " in " << Units::display(time, "ns") << " time, removing";
		trapped_charges_count += charge_per_step;
		}

		LOG(DEBUG) << " Propagated " << charge_per_step << " to " << Units::display(final_position, {"mm", "um"})
		<< " in " << Units::display(time, "ns") << " time, gain " << gain
		<< ", final state: " << allpix::to_string(state);

		// Create a new propagated charge and add it to the list
		auto global_position = detector_->getGlobalPosition(final_position);
		PropagatedCharge propagated_charge(final_position,
		global_position,
		deposit.getType(),
		static_cast<unsigned int>(std::lround(charge_per_step * gain)),
		deposit.getLocalTime() + time,
		deposit.getGlobalTime() + time,
		state,
		&deposit);

		propagated_charges.push_back(std::move(propagated_charge));
		deposit.getGlobalTime(),
		0,
		propagated_charges,
		output_plot_points);

		// Update statistical information
		++step_count;
		propagated_charges_count += charge_per_step;
		total_time += charge_per_step * time;
		if(output_plots_) {
		drift_time_histo_->Fill(static_cast<double>(Units::convert(time, "ns")), charge_per_step);
		group_size_histo_->Fill(charge_per_step);
		}
		recombined_charges_count += recombined;
		trapped_charges_count += trapped;
		propagated_charges_count += propagated;
		step_count += steps;
		total_time += time;
		}
		}

		@@ -413,18 +419,45 @@ void GenericPropagationModule::run(Event* event) {
		* velocity at every point with help of the electric field map of the detector. An Runge-Kutta integration is applied in
		* multiple steps, adding a random diffusion to the propagating charge every step.
		*/
		std::tuple<ROOT::Math::XYZPoint, double, double, CarrierState>
		GenericPropagationModule::propagate(const ROOT::Math::XYZPoint& pos,
		std::tuple<unsigned int, unsigned int, unsigned int, unsigned int, long double>
		GenericPropagationModule::propagate(Event* event,
		const DepositedCharge& deposit,
		const ROOT::Math::XYZPoint& pos,
		const CarrierType& type,
		const double initial_time,
		RandomNumberGenerator& random_generator,
		LineGraph::OutputPlotPoints& output_plot_points,
		const unsigned int charge) const {
		const unsigned int charge,
		const double initial_time_local,
		const double initial_time_global,
		const unsigned int level,
		std::vector<PropagatedCharge>& propagated_charges,
		LineGraph::OutputPlotPoints& output_plot_points) const {

		if(level > max_multiplication_level_) {
		LOG(WARNING) << "Found impact ionization shower with level larger than " << max_multiplication_level_
		<< ", interrupting";
		return {};
		}

		// Create a runge kutta solver using the electric field as step function
		Eigen::Vector3d position(pos.x(), pos.y(), pos.z());

		unsigned int propagated_charges_count = 0;
		unsigned int recombined_charges_count = 0;
		unsigned int trapped_charges_count = 0;
		unsigned int steps = 0;
		long double total_time = 0;

		// Add point of deposition to the output plots if requested
		if(output_linegraphs_) {
		output_plot_points.emplace_back(
		std::make_tuple(deposit.getGlobalTime(), charge, deposit.getType(), CarrierState::MOTION),
		std::vector<ROOT::Math::XYZPoint>());
		}
		auto output_plot_index = output_plot_points.size() - 1;

		// Initialize gain
		double gain = 1.;
		double gain_previous = 1.;
		unsigned int gain_integer = 1;

		// Define a function to compute the diffusion
		auto carrier_diffusion = [&](double efield_mag, double doping_concentration, double timestep) -> Eigen::Vector3d {
		@@ -433,9 +466,9 @@ GenericPropagationModule::propagate(const ROOT::Math::XYZPoint& pos,

		// Compute the independent diffusion in three
		allpix::normal_distribution<double> gauss_distribution(0, diffusion_std_dev);
		auto x = gauss_distribution(random_generator);
		auto y = gauss_distribution(random_generator);
		auto z = gauss_distribution(random_generator);
		auto x = gauss_distribution(event->getRandomEngine());
		auto y = gauss_distribution(event->getRandomEngine());
		auto z = gauss_distribution(event->getRandomEngine());
		return Eigen::Vector3d(x, y, z);
		};

		@@ -487,13 +520,13 @@ GenericPropagationModule::propagate(const ROOT::Math::XYZPoint& pos,
		double last_time = 0;
		size_t next_idx = 0;
		auto state = CarrierState::MOTION;
		while(state == CarrierState::MOTION && (initial_time + runge_kutta.getTime()) < integration_time_) {
		while(state == CarrierState::MOTION && (initial_time_local + runge_kutta.getTime()) < integration_time_) {
		// Update output plots if necessary (depending on the plot step)
		if(output_linegraphs_) {
		auto time_idx = static_cast<size_t>(runge_kutta.getTime() / output_plots_step_);
		while(next_idx <= time_idx) {
		output_plot_points.back().second.push_back(static_cast<ROOT::Math::XYZPoint>(position));
		next_idx = output_plot_points.back().second.size();
		output_plot_points.at(output_plot_index).second.push_back(static_cast<ROOT::Math::XYZPoint>(position));
		next_idx = output_plot_points.at(output_plot_index).second.size();
		}
		}

		@@ -529,19 +562,19 @@ GenericPropagationModule::propagate(const ROOT::Math::XYZPoint& pos,
		// Check if charge carrier is still alive:
		if(recombination_(type,
		detector_->getDopingConcentration(static_cast<ROOT::Math::XYZPoint>(position)),
		uniform_distribution(random_generator),
		uniform_distribution(event->getRandomEngine()),
		timestep)) {
		state = CarrierState::RECOMBINED;
		}

		// Check if the charge carrier has been trapped:
		if(trapping_(type, uniform_distribution(random_generator), timestep, std::sqrt(efield.Mag2()))) {
		if(trapping_(type, uniform_distribution(event->getRandomEngine()), timestep, std::sqrt(efield.Mag2()))) {
		if(output_plots_) {
		trapping_time_histo_->Fill(static_cast<double>(Units::convert(runge_kutta.getTime(), "ns")), charge);
		}

		auto detrap_time = detrapping_(type, uniform_distribution(random_generator), std::sqrt(efield.Mag2()));
		if((initial_time + runge_kutta.getTime() + detrap_time) < integration_time_) {
		auto detrap_time = detrapping_(type, uniform_distribution(event->getRandomEngine()), std::sqrt(efield.Mag2()));
		if((initial_time_local + runge_kutta.getTime() + detrap_time) < integration_time_) {
		LOG(DEBUG) << "De-trapping charge carrier after " << Units::display(detrap_time, {"ns", "us"});
		// De-trap and advance in time if still below integration time
		runge_kutta.advanceTime(detrap_time);
		@@ -557,18 +590,67 @@ GenericPropagationModule::propagate(const ROOT::Math::XYZPoint& pos,

		LOG(TRACE) << "Step from " << Units::display(static_cast<ROOT::Math::XYZPoint>(last_position), {"um", "mm"})
		<< " to " << Units::display(static_cast<ROOT::Math::XYZPoint>(position), {"um", "mm"}) << " at "
		<< Units::display(initial_time + runge_kutta.getTime(), {"ps", "ns", "us"})
		<< Units::display(initial_time_local + runge_kutta.getTime(), {"ps", "ns", "us"})
		<< ", state: " << allpix::to_string(state);

		// Apply multiplication step, fully deterministic from local efield and step length; Interpolate efield values
		gain *= multiplication_(type, (std::sqrt(efield.Mag2()) + std::sqrt(last_efield.Mag2())) / 2., step.value.norm());
		if(gain > 20.) {
		LOG(WARNING) << "Detected gain of " << gain << ", local electric field of "
		<< Units::display(std::sqrt(efield.Mag2()), "kV/cm") << ", diode seems to be in breakdown";
		} else if(gain > 1.) {
		if(gain > gain_previous) {
		LOG(DEBUG) << "Calculated gain of " << gain << " for step of " << Units::display(step.value.norm(), {"um", "nm"})
		<< " from field of " << Units::display(std::sqrt(last_efield.Mag2()), "kV/cm") << " to "
		<< Units::display(std::sqrt(efield.Mag2()), "kV/cm");

		// Update already processed gain factor:
		gain_previous = gain;

		if(gain > 20.) {
		LOG(WARNING) << "Detected gain of " << gain << ", local electric field of "
		<< Units::display(std::sqrt(efield.Mag2()), "kV/cm") << ", diode seems to be in breakdown";
		}

		// If the gain increased, we need to generate new charge carriers of the opposite type
		// Same-type carriers are simulated via the gain factor:
		auto floor_gain = static_cast<unsigned int>(std::floor(gain));

		auto inverted_type = invertCarrierType(type);
		auto do_propagate = [&]() {
		return (type == CarrierType::ELECTRON && propagate_electrons_) \|\|
		(type == CarrierType::HOLE && propagate_holes_);
		};

		if(gain_integer < floor_gain && do_propagate()) {
		// Placing new charge carrier mid-step:
		auto carrier_pos = static_cast<ROOT::Math::XYZPoint>(last_position + position) / 2.;
		LOG(DEBUG) << "Set of charge carriers (" << inverted_type << ") from gain on "
		<< Units::display(carrier_pos, {"mm", "um"});
		if(output_plots_) {
		multiplication_depth_histo_->Fill(carrier_pos.z(), charge * (floor_gain - gain_integer));
		}

		auto [recombined, trapped, propagated, psteps, ptime] =
		propagate(event,
		deposit,
		carrier_pos,
		inverted_type,
		charge * (floor_gain - gain_integer),
		initial_time_local + runge_kutta.getTime(),
		initial_time_global + runge_kutta.getTime(),
		level + 1,
		propagated_charges,
		output_plot_points);

		// Update statistics:
		recombined_charges_count += recombined;
		trapped_charges_count += trapped;
		propagated_charges_count += propagated;
		steps += psteps;
		total_time += ptime * charge;

		// Update the gain factor we have already generated opposite type charges for:
		gain_integer = floor_gain;
		LOG(DEBUG) << "Continuing propagation of charge carrier set (" << type << ") at "
		<< Units::display(carrier_pos, {"mm", "um"});
		}
		}

		// Update step length histogram
		@@ -612,14 +694,24 @@ GenericPropagationModule::propagate(const ROOT::Math::XYZPoint& pos,
		// If drift time is larger than integration time or the charge carriers have been collected at the backside, reset:
		if(!model_->isWithinImplant(static_cast<ROOT::Math::XYZPoint>(position)) &&
		(time >= integration_time_ \|\| last_position.z() < -model_->getSensorSize().z() * 0.45)) {
		std::get<3>(output_plot_points.back().first) = CarrierState::UNKNOWN;
		std::get<3>(output_plot_points.at(output_plot_index).first) = CarrierState::UNKNOWN;
		} else {
		std::get<3>(output_plot_points.back().first) = state;
		std::get<3>(output_plot_points.at(output_plot_index).first) = state;
		}
		}

		if(output_plots_) {
		gain_histo_->Fill(gain);
		if(output_plots_ && !multiplication_.is<NoImpactIonization>()) {
		if(level == 0) {
		gain_primary_histo_->Fill(gain, charge);
		if(type == CarrierType::ELECTRON) {
		gain_e_histo_->Fill(gain, charge);
		} else {
		gain_h_histo_->Fill(gain, charge);
		}
		}
		gain_all_histo_->Fill(gain, charge);

		multiplication_level_histo_->Fill(level, charge);
		}

		if(state == CarrierState::RECOMBINED) {
		@@ -629,8 +721,48 @@ GenericPropagationModule::propagate(const ROOT::Math::XYZPoint& pos,
		<< Units::display(static_cast<ROOT::Math::XYZPoint>(position), {"um", "mm"});
		}

		// Return the final position of the propagated charge, the time it took to propagate and its final state
		return std::make_tuple(static_cast<ROOT::Math::XYZPoint>(position), time, gain, state);
		auto local_position = static_cast<ROOT::Math::XYZPoint>(position);
		auto final_charge = static_cast<unsigned int>(charge * gain);

		if(state == CarrierState::RECOMBINED) {
		LOG(DEBUG) << " Recombined " << final_charge << " at " << Units::display(local_position, {"mm", "um"}) << " in "
		<< Units::display(time, "ns") << " time, removing";
		recombined_charges_count += final_charge;
		if(output_plots_) {
		recombination_time_histo_->Fill(static_cast<double>(Units::convert(time, "ns")), final_charge);
		}
		} else if(state == CarrierState::TRAPPED) {
		LOG(DEBUG) << " Trapped " << final_charge << " at " << Units::display(local_position, {"mm", "um"}) << " in "
		<< Units::display(time, "ns") << " time, removing";
		trapped_charges_count += final_charge;
		}
		propagated_charges_count += charge;
		++steps;
		total_time += time * charge;

		LOG(DEBUG) << " Propagated " << final_charge << " to " << Units::display(local_position, {"mm", "um"}) << " in "
		<< Units::display(time, "ns") << " time, gain " << gain << ", final state: " << allpix::to_string(state);

		// Create a new propagated charge and add it to the list
		auto global_position = detector_->getGlobalPosition(local_position);
		PropagatedCharge propagated_charge(local_position,
		global_position,
		deposit.getType(),
		static_cast<unsigned int>(std::lround(charge * gain)),
		deposit.getLocalTime() + time,
		deposit.getGlobalTime() + time,
		state,
		&deposit);

		propagated_charges.push_back(std::move(propagated_charge));

		if(output_plots_) {
		drift_time_histo_->Fill(static_cast<double>(Units::convert(time, "ns")), charge);
		group_size_histo_->Fill(charge);
		}

		// Return statistics counters about this and all daughter propagated charge carrier groups and their final states
		return std::make_tuple(recombined_charges_count, trapped_charges_count, propagated_charges_count, steps, total_time);
		}

		void GenericPropagationModule::finalize() {
		@@ -646,7 +778,14 @@ void GenericPropagationModule::finalize() {
		recombination_time_histo_->Write();
		trapping_time_histo_->Write();
		detrapping_time_histo_->Write();
		gain_histo_->Write();
		if(!multiplication_.is<NoImpactIonization>()) {
		gain_primary_histo_->Write();
		gain_all_histo_->Write();
		gain_e_histo_->Write();
		gain_h_histo_->Write();
		multiplication_level_histo_->Write();
		multiplication_depth_histo_->Write();
		}
		}

		long double average_time = static_cast<long double>(total_time_picoseconds_) / 1e3 /

src/modules/GenericPropagation/GenericPropagationModule.hpp

+28 −14

File changed.

Preview size limit exceeded, changes collapsed.

src/modules/GenericPropagation/README.md

+1 −0

Original line number	Diff line number	Diff line
		@@ -66,6 +66,7 @@ This module requires an installation of Eigen3.
		* `ignore_magnetic_field`: The magnetic field, if present, is ignored for this module. Defaults to false.
		* `multiplication_model`: Model used to calculate impact ionization parameters and charge multiplication. Defaults to `none` which corresponds to unity gain, a list of available models can be found in the documentation.
		* `multiplication_threshold`: Threshold field above which charge multiplication is calculated. Defaults to `100kV/cm`.
		* `max_multiplication_level`: Maximum level depth of the generated impact ionization charge multiplication shower after which the generation of further multiplication charge carrier levels is prohibited. This number represents the maximum number of daughter charge carrier groups that can be produced by one initial charge carrier group. This does not concern the size of the charge group itself but solely the level of generation. If a group generates a secondary group through impact ionization, the depth is `1`. If this secondary group again creates charge carriers when propagating, the level is `2` and so on. The default value is `5`.

		## Plotting parameters
		* `output_plots` : Determines if simple output plots should be generated for a monitoring of the simulation flow. Disabled by default.

src/modules/GenericPropagation/tests/13-impact-ionization.conf

+3 −0

Original line number	Diff line number	Diff line
		@@ -27,4 +27,7 @@ timestep_max = 1ps
		multiplication_model = "okuto"
		multiplication_threshold = 100kV/cm

		propagate_electrons = true
		propagate_holes = true

		#PASS Calculated gain of 1.10271 for step of 106.677nm from field of 186.529kV/cm to 186.646kV/cm

src/modules/GenericPropagation/tests/14-impact-missing-holes.conf

0 → 100644

+33 −0

Original line number	Diff line number	Diff line
		# SPDX-FileCopyrightText: 2022-2023 CERN and the Allpix Squared authors
		# SPDX-License-Identifier: MIT

		#DESC checks that a warning is printed in cases not all impact ionization carriers would be propagated
		[Allpix]
		detectors_file = "detector.conf"
		number_of_events = 1
		random_seed = 0

		[DepositionPointCharge]
		model = "fixed"
		source_type = "point"
		position = 445um 220um 50um
		number_of_charges = 1

		[ElectricFieldReader]
		model = "linear"
		bias_voltage = -1.4kV
		depletion_depth = 150um

		[GenericPropagation]
		log_level = DEBUG
		temperature = 293K
		charge_per_step = 1

		timestep_max = 1ps
		multiplication_model = "okuto"
		multiplication_threshold = 100kV/cm

		propagate_electrons = true
		propagate_holes = false

		#PASS (WARNING) [I:GenericPropagation:mydetector] Not propagating both types of charge carriers with charge multiplication enabled may lead to unphysical results