|
| 1 | +/** |
| 2 | + * @file |
| 3 | + * @brief Implementation of ProjectionPropagation module |
| 4 | + * @copyright MIT License |
| 5 | + */ |
| 6 | + |
| 7 | +#include "ProjectionPropagationModule.hpp" |
| 8 | + |
| 9 | +#include <cmath> |
| 10 | +#include <limits> |
| 11 | +#include <string> |
| 12 | +#include <utility> |
| 13 | + |
| 14 | +#include "core/messenger/Messenger.hpp" |
| 15 | +#include "core/utils/log.h" |
| 16 | +#include "objects/DepositedCharge.hpp" |
| 17 | +#include "objects/PropagatedCharge.hpp" |
| 18 | + |
| 19 | +using namespace allpix; |
| 20 | + |
| 21 | +ProjectionPropagationModule::ProjectionPropagationModule(Configuration config, |
| 22 | + Messenger* messenger, |
| 23 | + std::shared_ptr<Detector> detector) |
| 24 | + : Module(config, detector), config_(std::move(config)), messenger_(messenger), detector_(std::move(detector)) { |
| 25 | + // Save detector model |
| 26 | + model_ = detector_->getModel(); |
| 27 | + |
| 28 | + random_generator_.seed(getRandomSeed()); |
| 29 | + |
| 30 | + // Require deposits message for single detector |
| 31 | + messenger_->bindSingle(this, &ProjectionPropagationModule::deposits_message_, MsgFlags::REQUIRED); |
| 32 | + |
| 33 | + // Set default value for config variables |
| 34 | + config_.setDefault<int>("charge_per_step", 10); |
| 35 | + config_.setDefault<bool>("output_plots", false); |
| 36 | + |
| 37 | + output_plots_ = config_.get<bool>("output_plots"); |
| 38 | + |
| 39 | + // Set default for charge carrier propagation: |
| 40 | + config_.setDefault<bool>("propagate_holes", false); |
| 41 | + if(config_.get<bool>("propagate_holes")) { |
| 42 | + propagate_type_ = CarrierType::HOLE; |
| 43 | + LOG(INFO) << "Holes are chosen for propagation. Electrons are therefore not propagated."; |
| 44 | + } else { |
| 45 | + propagate_type_ = CarrierType::ELECTRON; |
| 46 | + } |
| 47 | + |
| 48 | + // Parameterization variables from https://doi.org/10.1016/0038-1101(77)90054-5 (section 5.2) |
| 49 | + auto temperature = config_.get<double>("temperature"); |
| 50 | + electron_Vm_ = Units::get(1.53e9 * std::pow(temperature, -0.87), "cm/s"); |
| 51 | + electron_Ec_ = Units::get(1.01 * std::pow(temperature, 1.55), "V/cm"); |
| 52 | + electron_Beta_ = 2.57e-2 * std::pow(temperature, 0.66); |
| 53 | + |
| 54 | + hole_Vm_ = Units::get(1.62e8 * std::pow(temperature, -0.52), "cm/s"); |
| 55 | + hole_Ec_ = Units::get(1.24 * std::pow(temperature, 1.68), "V/cm"); |
| 56 | + hole_Beta_ = 0.46 * std::pow(temperature, 0.17); |
| 57 | + |
| 58 | + boltzmann_kT_ = Units::get(8.6173e-5, "eV/K") * temperature; |
| 59 | +} |
| 60 | + |
| 61 | +void ProjectionPropagationModule::init() { |
| 62 | + if(detector_->getElectricFieldType() != ElectricFieldType::LINEAR) { |
| 63 | + throw ModuleError("This module should only be used with linear electric fields."); |
| 64 | + } |
| 65 | + |
| 66 | + if(output_plots_) { |
| 67 | + // Initialize output plot |
| 68 | + drift_time_histo = new TH1D("drift_time_histo", "Drift time;t[ns];particles", 75, 0., 25.); |
| 69 | + } |
| 70 | +} |
| 71 | + |
| 72 | +void ProjectionPropagationModule::run(unsigned int) { |
| 73 | + |
| 74 | + // Create vector of propagated charges to output |
| 75 | + std::vector<PropagatedCharge> propagated_charges; |
| 76 | + auto model = detector_->getModel(); |
| 77 | + |
| 78 | + double charge_lost = 0; |
| 79 | + double total_charge = 0; |
| 80 | + double total_projected_charge = 0; |
| 81 | + |
| 82 | + // Loop over all deposits for propagation |
| 83 | + for(auto& deposit : deposits_message_->getData()) { |
| 84 | + |
| 85 | + auto position = deposit.getLocalPosition(); |
| 86 | + auto type = deposit.getType(); |
| 87 | + |
| 88 | + // Selection of charge carrier: |
| 89 | + if(type != propagate_type_) { |
| 90 | + continue; |
| 91 | + } |
| 92 | + |
| 93 | + LOG(DEBUG) << "Set of " << deposit.getCharge() << " charge carriers (" << type << ") on " |
| 94 | + << display_vector(position, {"mm", "um"}); |
| 95 | + |
| 96 | + // Define a lambda function to compute the carrier mobility |
| 97 | + auto carrier_mobility = [&](double efield_mag) { |
| 98 | + // Compute carrier mobility from constants and electric field magnitude |
| 99 | + double numerator, denominator; |
| 100 | + if(type == CarrierType::ELECTRON) { |
| 101 | + numerator = electron_Vm_ / electron_Ec_; |
| 102 | + denominator = std::pow(1. + std::pow(efield_mag / electron_Ec_, electron_Beta_), 1.0 / electron_Beta_); |
| 103 | + } else { |
| 104 | + numerator = hole_Vm_ / hole_Ec_; |
| 105 | + denominator = std::pow(1. + std::pow(efield_mag / hole_Ec_, hole_Beta_), 1.0 / hole_Beta_); |
| 106 | + } |
| 107 | + return numerator / denominator; |
| 108 | + }; |
| 109 | + |
| 110 | + // Get the electric field at the position of the deposited charge and the top of the sensor: |
| 111 | + auto efield = detector_->getElectricField(position); |
| 112 | + double efield_mag = std::sqrt(efield.Mag2()); |
| 113 | + auto efield_top = detector_->getElectricField(ROOT::Math::XYZPoint(0., 0., model->getSensorSize().z() / 2.)); |
| 114 | + double efield_mag_top = std::sqrt(efield_top.Mag2()); |
| 115 | + |
| 116 | + LOG(TRACE) << "Electric field at carrier position / top of the sensor: " << Units::display(efield_mag_top, "V/cm") |
| 117 | + << " , " << Units::display(efield_mag, "V/cm"); |
| 118 | + |
| 119 | + slope_efield_ = (efield_mag_top - efield_mag) / (model->getSensorSize().z() / 2. - position.z()); |
| 120 | + |
| 121 | + // Calculate the drift time |
| 122 | + auto calc_drift_time = [&]() { |
| 123 | + double Ec = (type == CarrierType::ELECTRON ? electron_Ec_ : hole_Ec_); |
| 124 | + double zero_mobility = (type == CarrierType::ELECTRON ? electron_Vm_ / electron_Ec_ : hole_Vm_ / hole_Ec_); |
| 125 | + |
| 126 | + return ((log(efield_mag_top) - log(efield_mag)) / slope_efield_ + |
| 127 | + (model->getSensorSize().z() / 2. - position.z()) / Ec) / |
| 128 | + zero_mobility; |
| 129 | + }; |
| 130 | + |
| 131 | + // Only project if within the depleted region (i.e. efield not zero) |
| 132 | + if(efield_mag < std::numeric_limits<double>::epsilon()) { |
| 133 | + LOG(TRACE) << "Electric field is zero at " << display_vector(position, {"mm", "um"}); |
| 134 | + continue; |
| 135 | + } |
| 136 | + |
| 137 | + LOG(TRACE) << "Electric field is " << Units::display(efield_mag, "V/cm"); |
| 138 | + |
| 139 | + // Assume linear electric field over the sensor: |
| 140 | + double diffusion_constant = boltzmann_kT_ * (carrier_mobility(efield_mag) + carrier_mobility(efield_mag_top)) / 2.; |
| 141 | + |
| 142 | + double drift_time = calc_drift_time(); |
| 143 | + LOG(TRACE) << "Drift time is " << Units::display(drift_time, "ns"); |
| 144 | + |
| 145 | + if(output_plots_) { |
| 146 | + drift_time_histo->SetBinContent(drift_time_histo->FindBin(drift_time), |
| 147 | + drift_time_histo->GetBinContent(drift_time_histo->FindBin(drift_time)) + |
| 148 | + deposit.getCharge()); |
| 149 | + } |
| 150 | + |
| 151 | + double diffusion_std_dev = std::sqrt(2. * diffusion_constant * drift_time); |
| 152 | + LOG(TRACE) << "Diffusion width is " << Units::display(diffusion_std_dev, "um"); |
| 153 | + |
| 154 | + double projected_charge = 0; |
| 155 | + |
| 156 | + unsigned int charges_remaining = deposit.getCharge(); |
| 157 | + total_charge += charges_remaining; |
| 158 | + |
| 159 | + auto charge_per_step = config_.get<unsigned int>("charge_per_step"); |
| 160 | + while(charges_remaining > 0) { |
| 161 | + if(charge_per_step > charges_remaining) { |
| 162 | + charge_per_step = charges_remaining; |
| 163 | + } |
| 164 | + charges_remaining -= charge_per_step; |
| 165 | + |
| 166 | + std::normal_distribution<double> gauss_distribution(0, diffusion_std_dev); |
| 167 | + double diffusion_x = gauss_distribution(random_generator_); |
| 168 | + double diffusion_y = gauss_distribution(random_generator_); |
| 169 | + |
| 170 | + auto projected_position = ROOT::Math::XYZPoint( |
| 171 | + position.x() + diffusion_x, position.y() + diffusion_y, model->getSensorSize().z() / 2.); |
| 172 | + |
| 173 | + // Only add if within sensor volume: |
| 174 | + auto local_position = |
| 175 | + ROOT::Math::XYZPoint(projected_position.x(), projected_position.y(), model->getSensorSize().z() / 2.); |
| 176 | + if(!detector_->isWithinSensor(local_position)) { |
| 177 | + continue; |
| 178 | + } |
| 179 | + |
| 180 | + auto global_position = detector_->getGlobalPosition(local_position); |
| 181 | + |
| 182 | + // Produce charge carrier at this position |
| 183 | + propagated_charges.emplace_back( |
| 184 | + projected_position, global_position, deposit.getType(), charge_per_step, deposit.getEventTime(), &deposit); |
| 185 | + |
| 186 | + LOG(DEBUG) << "Propagated " << charge_per_step << " charge carriers (" << type << ") to " |
| 187 | + << display_vector(projected_position, {"mm", "um"}); |
| 188 | + |
| 189 | + projected_charge += charge_per_step; |
| 190 | + } |
| 191 | + total_projected_charge += projected_charge; |
| 192 | + } |
| 193 | + charge_lost = total_charge - total_projected_charge; |
| 194 | + |
| 195 | + LOG(INFO) << "Total charge: " << total_charge << " (lost: " << charge_lost << ", " << (charge_lost / total_charge * 100.) |
| 196 | + << "%)"; |
| 197 | + LOG(DEBUG) << "Total count of propagated charge carriers: " << propagated_charges.size(); |
| 198 | + |
| 199 | + // Create a new message with propagated charges |
| 200 | + auto propagated_charge_message = std::make_shared<PropagatedChargeMessage>(std::move(propagated_charges), detector_); |
| 201 | + |
| 202 | + // Dispatch the message with propagated charges |
| 203 | + messenger_->dispatchMessage(this, propagated_charge_message); |
| 204 | +} |
| 205 | + |
| 206 | +void ProjectionPropagationModule::finalize() { |
| 207 | + if(output_plots_) { |
| 208 | + // Write output plot |
| 209 | + drift_time_histo->Write(); |
| 210 | + } |
| 211 | +} |
0 commit comments