Commit 77e42a40 authored by Simon Spannagel's avatar Simon Spannagel
Browse files

GenericPropagation: use CarrierState to govern propagation loop

parent c9c9de7c
Loading
Loading
Loading
Loading
+21 −15
Original line number Diff line number Diff line
@@ -623,10 +623,10 @@ void GenericPropagationModule::run(Event* event) {
            }

            // Propagate a single charge deposit
            auto [final_position, time, alive] = propagate(
            auto [final_position, time, state] = propagate(
                initial_position, deposit.getType(), deposit.getLocalTime(), event->getRandomEngine(), output_plot_points);

            if(!alive) {
            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;
@@ -644,7 +644,7 @@ void GenericPropagationModule::run(Event* event) {
                                               charge_per_step,
                                               deposit.getLocalTime() + time,
                                               deposit.getGlobalTime() + time,
                                               (alive ? CarrierState::MOTION : CarrierState::RECOMBINED),
                                               state,
                                               &deposit);

            propagated_charges.push_back(std::move(propagated_charge));
@@ -691,7 +691,7 @@ 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, bool>
std::tuple<ROOT::Math::XYZPoint, double, CarrierState>
GenericPropagationModule::propagate(const ROOT::Math::XYZPoint& pos,
                                    const CarrierType& type,
                                    const double initial_time,
@@ -759,9 +759,8 @@ GenericPropagationModule::propagate(const ROOT::Math::XYZPoint& pos,
    Eigen::Vector3d last_position = position;
    double last_time = 0;
    size_t next_idx = 0;
    bool is_alive = true;
    while(detector_->getModel()->isWithinSensor(static_cast<ROOT::Math::XYZPoint>(position)) &&
          (initial_time + runge_kutta.getTime()) < integration_time_ && is_alive) {
    auto state = CarrierState::MOTION;
    while(state == CarrierState::MOTION && (initial_time + 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_);
@@ -791,16 +790,23 @@ GenericPropagationModule::propagate(const ROOT::Math::XYZPoint& pos,
        position += diffusion;
        runge_kutta.setValue(position);

        // Check if we are still in the sensor:
        if(!detector_->getModel()->isWithinSensor(static_cast<ROOT::Math::XYZPoint>(position))) {
            state = CarrierState::HALTED;
        }

        // Check if charge carrier is still alive:
        is_alive = !recombination_(type,
        if(recombination_(type,
                          detector_->getDopingConcentration(static_cast<ROOT::Math::XYZPoint>(position)),
                          survival(random_generator),
                                   timestep);
                          timestep)) {
            state = CarrierState::RECOMBINED;
        }

        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"})
                   << (is_alive ? "" : ", recombined");
                   << ", state: " << allpix::to_string(state);
        // Adapt step size to match target precision
        double uncertainty = step.error.norm();

@@ -831,7 +837,7 @@ GenericPropagationModule::propagate(const ROOT::Math::XYZPoint& pos,

    // Find proper final position in the sensor
    auto time = runge_kutta.getTime();
    if(!detector_->getModel()->isWithinSensor(static_cast<ROOT::Math::XYZPoint>(position))) {
    if(state == CarrierState::HALTED) {
        auto check_position = position;
        check_position.z() = last_position.z();
        if(position.z() > 0 && detector_->getModel()->isWithinSensor(static_cast<ROOT::Math::XYZPoint>(check_position))) {
@@ -856,12 +862,12 @@ GenericPropagationModule::propagate(const ROOT::Math::XYZPoint& pos,
        }
    }

    if(!is_alive) {
    if(state == CarrierState::RECOMBINED) {
        LOG(DEBUG) << "Charge carrier recombined after " << Units::display(last_time, {"ns"});
    }

    // Return the final position of the propagated charge
    return std::make_tuple(static_cast<ROOT::Math::XYZPoint>(position), initial_time + time, is_alive);
    return std::make_tuple(static_cast<ROOT::Math::XYZPoint>(position), initial_time + time, state);
}

void GenericPropagationModule::finalize() {
+7 −7
Original line number Diff line number Diff line
@@ -88,10 +88,10 @@ namespace allpix {
         * @param initial_time Initial time passed before propagation starts in local time coordinates
         * @param random_generator Reference to the random number engine to be used
         * @param output_plot_points Reference to vector to hold points for line graph output plots
         * @return Tuple with the point where the deposit ended after propagation, the time the propagation took and a flag
         * whether it has recombined
         * @return Tuple with the point where the deposit ended after propagation, the time the propagation took and the
         * final state of the charge carrier at the end of processing
         */
        std::tuple<ROOT::Math::XYZPoint, double, bool> propagate(const ROOT::Math::XYZPoint& pos,
        std::tuple<ROOT::Math::XYZPoint, double, CarrierState> propagate(const ROOT::Math::XYZPoint& pos,
                                                                         const CarrierType& type,
                                                                         const double initial_time,
                                                                         RandomNumberGenerator& random_generator,