ImpactIonisation: apply changes made to TransProp also to GenProp (09b08a9b) · Commits · XRFLab109 / Allpix Squared

src/modules/GenericPropagation/GenericPropagationModule.cpp

+55 −46

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,
		@@ -455,9 +455,7 @@ 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;
		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 transverse 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 from a geometric distribution (Yule process) 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::floor(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,