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

doc/usermanual/06_models/05_impact_ionization.md

+18 −1

Original line number	Diff line number	Diff line
		@@ -8,7 +8,7 @@ weight: 5
		Allpix Squared implements charge multiplication via impact ionization models. These models are only used by propagation
		modules which perform a step-by-step simulation of the charge carrier motion.

		The gain $`g`$ is calculated for all models as exponential of the model-dependent impact ionization coefficient $`\alpha`$ and
		The per-step gain $`g`$ is calculated for all models as exponential of the model-dependent impact ionization coefficient $`\alpha`$ and
		the length of the step $`l`$ performed in the respective electric field. If the electric field strength stays below a
		configurable threshold $`E_{\text{thr}}`$, unity gain is assumed:

		@@ -21,6 +21,22 @@ g (E, T) = \left\{
		\right.
		```

		The impact ionization coefficient $`\alpha`$ is calculated depending on the selected impact ionization model. The models themselves are described below.

		The number of additional charge carriers generated per step $`n`$ is determined via a stochastic approach by applying the following equation dependent on a random number drawn from a uniform distribution $`u(0,1)`$
		```math
		n = \frac{\ln(u)}{\ln(1-1/g)} = \frac{1}{\log_u(1-1/g)}
		```
		This distribution is applied e.g. in Garfield++\[[@garfieldpp]\] and represents a microscopic simulation of Yule processes.

		The number of secondary charge carriers generated from impact ionization is calculated for every individual charge carrier within a group of charge carriers and summed per propagation step. Additional charge carriers are then added to the group (same-type carriers) or deposited (opposite-type) at the end of the corresponding step.

		This algorithm results in a mean number of secondaries generated equal to
		```math
		<n_{total}> = \exp\left(\int_{x_0}^{x_n}\alpha(x)dx \right)
		```
		for sufficiently low step sizes.

		The following impact ionization models are available:

		## Massey Model
		@@ -305,6 +321,7 @@ The interpretation of the custom impact ionization functions is based on the `RO
		supports all corresponding features, mathematical expressions and constants.


		[@garfieldpp]: https://gitlab.cern.ch/garfield/garfieldpp
		[@massey]: https://doi.org/10.1109/TED.2006.881010
		[@rd50ionization]: https://arxiv.org/abs/2211.16543
		[@overstraeten]: https://doi.org/10.1016/0038-1101(70)90139-5

src/modules/GenericPropagation/GenericPropagationModule.cpp

+56 −47

Original line number	Diff line number	Diff line
		@@ -233,25 +233,25 @@ void GenericPropagationModule::initialize() {
		gain_primary_histo_ = CreateHistogram<TH1D>(
		"gain_primary_histo",
		"Gain per primarily induced charge carrier group after propagation;gain;number of groups transported",
		500,
		24,
		1,
		25);
		gain_all_histo_ =
		CreateHistogram<TH1D>("gain_all_histo",
		"Gain per charge carrier group after propagation;gain;number of groups transported",
		500,
		24,
		1,
		25);
		gain_e_histo_ =
		CreateHistogram<TH1D>("gain_e_histo",
		"Gain per primary electron group after propagation;gain;number of groups transported",
		500,
		24,
		1,
		25);
		gain_h_histo_ =
		CreateHistogram<TH1D>("gain_h_histo",
		"Gain per primary hole group after propagation;gain;number of groups transported",
		500,
		24,
		1,
		25);
		multiplication_level_histo_ = CreateHistogram<TH1D>(
		@@ -424,7 +424,7 @@ GenericPropagationModule::propagate(Event* event,
		const DepositedCharge& deposit,
		const ROOT::Math::XYZPoint& pos,
		const CarrierType& type,
		const unsigned int charge,
		unsigned int charge,
		const double initial_time_local,
		const double initial_time_global,
		const unsigned int level,
		@@ -454,10 +454,8 @@ GenericPropagationModule::propagate(Event* event,
		}
		auto output_plot_index = output_plot_points.size() - 1;

		// Initialize gain
		double gain = 1.;
		double gain_previous = 1.;
		unsigned int gain_integer = 1;
		// Store initial charge
		const unsigned int initial_charge = charge;

		// Define a function to compute the diffusion
		auto carrier_diffusion = [&](double efield_mag, double doping_concentration, double timestep) -> Eigen::Vector3d {
		@@ -593,38 +591,43 @@ GenericPropagationModule::propagate(Event* event,
		<< 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 > 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");
		// Apply multiplication step: calculate gain factor from local efield and step length; Interpolate efield values
		// The multiplication factor is not scaled by the velocity fraction parallel to the electric field, as the
		// correction is negligible for semiconductors
		auto local_gain =
		multiplication_(type, (std::sqrt(efield.Mag2()) + std::sqrt(last_efield.Mag2())) / 2., step.value.norm());

		// Update already processed gain factor:
		gain_previous = gain;
		unsigned int n_secondaries = 0;

		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(local_gain > 1.0) {
		LOG(DEBUG) << "Calculated local gain of " << local_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");

		// 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));
		// For each charge carrier draw a number to determine the number of
		// secondaries generated in this step
		double log_prob = 1. / std::log1p(-1. / local_gain);
		for(unsigned int i_carrier = 0; i_carrier < charge; ++i_carrier) {
		n_secondaries +=
		static_cast<unsigned int>(std::log(uniform_distribution(event->getRandomEngine())) * log_prob);
		}

		auto inverted_type = invertCarrierType(type);
		auto do_propagate = [&]() {
		return (type == CarrierType::ELECTRON && propagate_electrons_) \|\|
		(type == CarrierType::HOLE && propagate_holes_);
		return (inverted_type == CarrierType::ELECTRON && propagate_electrons_) \|\|
		(inverted_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 "
		if(n_secondaries > 0 && do_propagate()) {
		// Generate new charge carriers of the opposite type
		// Same-type charge carriers are generated by increasing the charge at the end of the step
		// Placing new charge carrier at the end of the step:
		auto carrier_pos = static_cast<ROOT::Math::XYZPoint>(position);
		LOG(DEBUG) << "Set of charge carriers (" << inverted_type << ") generated from impact ionization on "
		<< Units::display(carrier_pos, {"mm", "um"});
		if(output_plots_) {
		multiplication_depth_histo_->Fill(carrier_pos.z(), charge * (floor_gain - gain_integer));
		multiplication_depth_histo_->Fill(carrier_pos.z(), n_secondaries);
		}

		auto [recombined, trapped, propagated, psteps, ptime] =
		@@ -632,7 +635,7 @@ GenericPropagationModule::propagate(Event* event,
		deposit,
		carrier_pos,
		inverted_type,
		charge * (floor_gain - gain_integer),
		n_secondaries,
		initial_time_local + runge_kutta.getTime(),
		initial_time_global + runge_kutta.getTime(),
		level + 1,
		@@ -646,11 +649,15 @@ GenericPropagationModule::propagate(Event* event,
		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"});
		}

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

		// Update step length histogram
		@@ -679,6 +686,8 @@ GenericPropagationModule::propagate(Event* event,
		timestep = timestep_min_;
		}
		runge_kutta.setTimeStep(timestep);

		charge += n_secondaries;
		}

		// Find proper final position in the sensor
		@@ -700,18 +709,19 @@ GenericPropagationModule::propagate(Event* event,
		}
		}

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

		multiplication_level_histo_->Fill(level, charge);
		multiplication_level_histo_->Fill(level, initial_charge);
		}

		if(state == CarrierState::RECOMBINED) {
		@@ -722,25 +732,24 @@ GenericPropagationModule::propagate(Event* event,
		}

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

		LOG(DEBUG) << " Propagated " << final_charge << " to " << Units::display(local_position, {"mm", "um"}) << " in "
		LOG(DEBUG) << " Propagated " << 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
		@@ -748,7 +757,7 @@ GenericPropagationModule::propagate(Event* event,
		PropagatedCharge propagated_charge(local_position,
		global_position,
		deposit.getType(),
		static_cast<unsigned int>(std::lround(charge * gain)),
		charge,
		deposit.getLocalTime() + time,
		deposit.getGlobalTime() + time,
		state,

src/modules/GenericPropagation/GenericPropagationModule.hpp

+1 −1

Original line number	Diff line number	Diff line
		@@ -98,7 +98,7 @@ namespace allpix {
		const DepositedCharge& deposit,
		const ROOT::Math::XYZPoint& pos,
		const CarrierType& type,
		const unsigned int charge,
		unsigned int charge,
		const double initial_time_local,
		const double initial_time_global,
		const unsigned int level,

src/modules/GenericPropagation/README.md

+2 −0

Original line number	Diff line number	Diff line
		@@ -13,7 +13,9 @@ module_output: "PropagatedCharge"
		Simulates the propagation of electrons and/or holes through the sensitive sensor volume of the detector. It allows to propagate sets of charge carriers together in order to speed up the simulation while maintaining the required accuracy. The propagation process for these sets is fully independent and no interaction is simulated. The maximum size of the set of propagated charges and thus the accuracy of the propagation can be controlled via the `charge_per_step` parameter. The maximum number of charge groups to be propagated for a single deposit position can be controlled via the `max_charge_groups` parameter.

		The propagation consists of a combination of drift and diffusion simulation. The drift is calculated using the charge carrier velocity derived from the charge carrier mobility and the magnetic field via a calculation of the Lorentz drift. The correct mobility for either electrons or holes is automatically chosen, based on the type of the charge carrier under consideration. Thus, also input with both electrons and holes is treated properly. The mobility model can be chosen using the `mobility_model` parameter, and a list of available models can be found in the user manual.

		This module implements charge multiplication by impact ionization. The multiplication model can be chosen using the `multiplication_model` parameter, the list of available models can be found in the user manual. By default, the model defaults to `none` and impact ionization is switched off, generating unity gain.
		To simulate impact ionization, the number of newly generated electron-hole pairs is calculated for every propagation step and every charge carrier in the group, based on drawing a random number from a geometric distribution. This represents a stepwise approach to the avalanche generation process. The charge of a charge group is increased by the number of impact ionization processes per step and opposite-type charge carriers are generated at the end of the step, if the opposite-type charge carrier is selected to be propagated (see below).

		The two parameters `propagate_electrons` and `propagate_holes` allow to control which type of charge carrier is propagated to their respective electrodes. Either one of the carrier types can be selected, or both can be propagated. It should be noted that this will slow down the simulation considerably since twice as many carriers have to be handled and it should only be used where sensible.
		The direction of the propagation depends on the electric and magnetic fields field configured, and it should be ensured that the carrier types selected are actually transported to the implant side. For linear electric fields, a warning is issued if a possible misconfiguration is detected.

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

+1 −1

Original line number	Diff line number	Diff line
		@@ -30,4 +30,4 @@ 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
		#PASS Calculated local gain of 1.00116 for step of 106.676nm from field of 186.385kV/cm to 186.547kV/cm