g4simulation/g4tpc/TpcPrimaryIonizationModel.cc

0001 #include "TpcPrimaryIonizationModel.h"
0002
0003 #include <algorithm>
0004 #include <array>
0005 #include <cmath>
0006 #include <iterator>
0007 #include <stdexcept>
0008 #include <vector>
0009
0010 namespace
0011 {
0012   constexpr double argon_primary_clusters_per_cm = 23.0;
0013   constexpr double cf4_primary_clusters_per_cm = 51.0;
0014   constexpr double isobutane_primary_clusters_per_cm = 84.0;
0015
0016   // Measured single-cluster electron multiplicities for Ar and CH4:
0017   // H. Fischle, J. Heintze, B. Schmidt, Nucl. Instrum. Meth. A 301
0018   // (1991) 202-214, doi:10.1016/0168-9002(91)90460-8.
0019   constexpr std::array<double, 19> argon_cluster_probabilities = {
0020       0.656, 0.150, 0.064, 0.035, 0.0225, 0.0155, 0.0105,
0021       0.0081, 0.0061, 0.0049, 0.0039, 0.0030, 0.0025, 0.0020,
0022       0.0016, 0.0012, 0.00095, 0.00075, 0.00063};
0023
0024   constexpr std::array<double, 19> methane_cluster_probabilities = {
0025       0.786, 0.120, 0.034, 0.016, 0.0095, 0.0060, 0.0044,
0026       0.0034, 0.0027, 0.0021, 0.0017, 0.0013, 0.0010, 0.0008,
0027       0.0006, 0.00050, 0.00042, 0.00037, 0.00033};
0028
0029   template <std::size_t Size>
0030   std::vector<double> build_cluster_probabilities(
0031       const std::array<double, Size> &measured)
0032   {
0033     std::vector<double> probabilities(
0034         TpcPrimaryIonizationModel::maximum_cluster_size + 1, 0.0);
0035
0036     double measured_sum = 0.0;
0037     for (std::size_t index = 0; index < measured.size(); ++index)
0038     {
0039       probabilities[index + 1] = measured[index];
0040       measured_sum += measured[index];
0041     }
0042
0043     if (!(measured_sum > 0.0 && measured_sum <= 1.0))
0044     {
0045       throw std::logic_error("invalid measured TPC cluster-size probabilities");
0046     }
0047
0048     // The published table is used through n=19. Allocate its remaining
0049     // probability to a truncated 1/n^2 tail through the explicit cutoff.
0050     constexpr unsigned int tail_begin = Size + 1;
0051     double inverse_square_sum = 0.0;
0052     for (unsigned int size = tail_begin;
0053          size <= TpcPrimaryIonizationModel::maximum_cluster_size;
0054          ++size)
0055     {
0056       inverse_square_sum += 1.0 / static_cast<double>(size * size);
0057     }
0058
0059     const double tail_normalization =
0060         (1.0 - measured_sum) / inverse_square_sum;
0061     for (unsigned int size = tail_begin;
0062          size <= TpcPrimaryIonizationModel::maximum_cluster_size;
0063          ++size)
0064     {
0065       probabilities[size] =
0066           tail_normalization / static_cast<double>(size * size);
0067     }
0068
0069     return probabilities;
0070   }
0071 }  // namespace
0072
0073 TpcPrimaryIonizationModel::TpcPrimaryIonizationModel(
0074     const GasFractions &fractions)
0075 {
0076   const std::array<double, 3> values = {
0077       fractions.argon, fractions.cf4, fractions.isobutane};
0078   for (const double fraction : values)
0079   {
0080     if (!std::isfinite(fraction) || fraction < 0.0)
0081     {
0082       throw std::invalid_argument(
0083           "TPC primary-ionization gas fractions must be finite and nonnegative");
0084     }
0085   }
0086
0087   const double fraction_sum =
0088       fractions.argon + fractions.cf4 + fractions.isobutane;
0089   if (!std::isfinite(fraction_sum) || std::abs(fraction_sum - 1.0) > 1.e-6)
0090   {
0091     throw std::invalid_argument(
0092         "TPC primary-ionization gas fractions must sum to one");
0093   }
0094
0095   const double argon_weight =
0096       fractions.argon * argon_primary_clusters_per_cm;
0097   const double cf4_weight = fractions.cf4 * cf4_primary_clusters_per_cm;
0098   const double isobutane_weight =
0099       fractions.isobutane * isobutane_primary_clusters_per_cm;
0100   m_primaryClustersPerCm = argon_weight + cf4_weight + isobutane_weight;
0101
0102   if (!(m_primaryClustersPerCm > 0.0))
0103   {
0104     throw std::invalid_argument(
0105         "TPC primary-ionization gas mixture has zero cluster density");
0106   }
0107
0108   const auto argon_probabilities =
0109       build_cluster_probabilities(argon_cluster_probabilities);
0110   const auto methane_probabilities =
0111       build_cluster_probabilities(methane_cluster_probabilities);
0112
0113   // No validated CF4 or isobutane table is supplied by this first model.
0114   // CF4 therefore uses the Ar shape and isobutane uses CH4 as an explicit
0115   // proxy. Each component is weighted by its contribution to the primary
0116   // cluster density rather than by volume fraction alone.
0117   const double argon_shape_weight =
0118       (argon_weight + cf4_weight) / m_primaryClustersPerCm;
0119   const double methane_shape_weight =
0120       isobutane_weight / m_primaryClustersPerCm;
0121
0122   m_clusterSizeCdf.resize(maximum_cluster_size + 1, 0.0);
0123   double cumulative_probability = 0.0;
0124   for (unsigned int size = 1; size <= maximum_cluster_size; ++size)
0125   {
0126     cumulative_probability +=
0127         argon_shape_weight * argon_probabilities[size] +
0128         methane_shape_weight * methane_probabilities[size];
0129     m_clusterSizeCdf[size] = cumulative_probability;
0130   }
0131   m_clusterSizeCdf[maximum_cluster_size] = 1.0;
0132 }
0133
0134 double TpcPrimaryIonizationModel::mean_primary_clusters(
0135     const double path_length_cm) const
0136 {
0137   if (!std::isfinite(path_length_cm) || path_length_cm < 0.0)
0138   {
0139     throw std::invalid_argument(
0140         "TPC primary-ionization path length must be finite and nonnegative");
0141   }
0142   return m_primaryClustersPerCm * path_length_cm;
0143 }
0144
0145 unsigned int TpcPrimaryIonizationModel::sample_cluster_size(
0146     const double uniform_01) const
0147 {
0148   if (!std::isfinite(uniform_01) || uniform_01 < 0.0 || uniform_01 >= 1.0)
0149   {
0150     throw std::invalid_argument(
0151         "TPC cluster-size random variate must lie in [0, 1)");
0152   }
0153
0154   const auto found = std::lower_bound(
0155       std::next(m_clusterSizeCdf.begin()),
0156       m_clusterSizeCdf.end(),
0157       uniform_01);
0158   return static_cast<unsigned int>(
0159       std::distance(m_clusterSizeCdf.begin(), found));
0160 }
0161
0162 double TpcPrimaryIonizationModel::mean_cluster_size() const
0163 {
0164   double previous_cdf = 0.0;
0165   double mean = 0.0;
0166   for (unsigned int size = 1; size <= maximum_cluster_size; ++size)
0167   {
0168     const double probability = m_clusterSizeCdf[size] - previous_cdf;
0169     mean += static_cast<double>(size) * probability;
0170     previous_cdf = m_clusterSizeCdf[size];
0171   }
0172   return mean;
0173 }