Interactions.cpp

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Interactions.cpp

// This file is part of the Acts project.
//
// Copyright (C) 2019-2020 CERN for the benefit of the Acts project
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.

#include "Acts/Material/Interactions.hpp"

#include "Acts/Definitions/PdgParticle.hpp"
#include "Acts/Material/Material.hpp"

#include <cassert>
#include <cmath>

using namespace Acts::UnitLiterals;

namespace {

// values from RPP2018 table 33.1
// electron mass
constexpr float Me = 0.5109989461_MeV;
// Bethe formular prefactor. 1/mol unit is just a factor 1 here.
constexpr float K = 0.307075_MeV * 1_cm * 1_cm;
// Energy scale for plasma energy.
constexpr float PlasmaEnergyScale = 28.816_eV;

struct RelativisticQuantities {
  float q2OverBeta2 = 0.0f;
  float beta2 = 0.0f;
  float betaGamma = 0.0f;
  float gamma = 0.0f;

  RelativisticQuantities(float mass, float qOverP, float absQ) {
    assert((0 < mass) and "Mass must be positive");
    assert((qOverP != 0) and "q/p must be non-zero");
    assert((absQ > 0) and "Absolute charge must be non-zero and positive");
    // beta²/q² = (p/E)²/q² = p²/(q²m² + q²p²) = 1/(q² + (m²(q/p)²)
    // q²/beta² = q² + m²(q/p)²
    q2OverBeta2 = absQ * absQ + (mass * qOverP) * (mass * qOverP);
    assert((q2OverBeta2 >= 0) && "Negative q2OverBeta2");
    // 1/p = q/(qp) = (q/p)/q
    const float mOverP = mass * std::abs(qOverP / absQ);
    const float pOverM = 1.0f / mOverP;
    // beta² = p²/E² = p²/(m² + p²) = 1/(1 + (m/p)²)
    beta2 = 1.0f / (1.0f + mOverP * mOverP);
    assert((beta2 >= 0) && "Negative beta2");
    // beta*gamma = (p/sqrt(m² + p²))*(sqrt(m² + p²)/m) = p/m
    betaGamma = pOverM;
    assert((betaGamma >= 0) && "Negative betaGamma");
    // gamma = sqrt(m² + p²)/m = sqrt(1 + (p/m)²)
    gamma = std::sqrt(1.0f + pOverM * pOverM);
  }
};

inline float deriveBeta2(float qOverP, const RelativisticQuantities& rq) {
  return -2 / (qOverP * rq.gamma * rq.gamma);
}

inline float computeMassTerm(float mass, const RelativisticQuantities& rq) {
  return 2 * mass * rq.betaGamma * rq.betaGamma;
}

inline float logDeriveMassTerm(float qOverP) {
  // only need to compute d((beta*gamma)²)/(beta*gamma)²; rest cancels.
  return -2 / qOverP;
}

inline float computeWMax(float mass, const RelativisticQuantities& rq) {
  const float mfrac = Me / mass;
  const float nominator = 2 * Me * rq.betaGamma * rq.betaGamma;
  const float denonimator = 1.0f + 2 * rq.gamma * mfrac + mfrac * mfrac;
  return nominator / denonimator;
}

inline float logDeriveWMax(float mass, float qOverP,
                           const RelativisticQuantities& rq) {
  // this is (q/p) * (beta/q).
  // both quantities have the same sign and the product must always be
  // positive. we can thus reuse the known (unsigned) quantity (q/beta)².
  const float a = std::abs(qOverP / std::sqrt(rq.q2OverBeta2));
  // (m² + me²) / me = me (1 + (m/me)²)
  const float b = Me * (1.0f + (mass / Me) * (mass / Me));
  return -2 * (a * b - 2 + rq.beta2) / (qOverP * (a * b + 2));
}

inline float computeEpsilon(float molarElectronDensity, float thickness,
                            const RelativisticQuantities& rq) {
  return 0.5f * K * molarElectronDensity * thickness * rq.q2OverBeta2;
}

inline float logDeriveEpsilon(float qOverP, const RelativisticQuantities& rq) {
  // only need to compute d(q²/beta²)/(q²/beta²); everything else cancels.
  return 2 / (qOverP * rq.gamma * rq.gamma);
}

inline float computeDeltaHalf(float meanExitationPotential,
                              float molarElectronDensity,
                              const RelativisticQuantities& rq) {
  // only relevant for very high ernergies; use arbitrary cutoff
  if (rq.betaGamma < 10.0f) {
    return 0.0f;
  }
  // pre-factor according to RPP2019 table 33.1
  const float plasmaEnergy =
      PlasmaEnergyScale * std::sqrt(1000.f * molarElectronDensity);
  return std::log(rq.betaGamma) +
         std::log(plasmaEnergy / meanExitationPotential) - 0.5f;
}

inline float deriveDeltaHalf(float qOverP, const RelativisticQuantities& rq) {
  // original equation is of the form
  //     log(beta*gamma) + log(eplasma/I) - 1/2
  // which the resulting derivative as
  //     d(beta*gamma)/(beta*gamma)
  return (rq.betaGamma < 10.0f) ? 0.0f : (-1.0f / qOverP);
}

inline float convertLandauFwhmToGaussianSigma(float fwhm) {
  // return fwhm / (2 * std::sqrt(2 * std::log(2.0f)));
  return fwhm * 0.42466090014400953f;
}

namespace detail {

inline float computeEnergyLossLandauFwhm(const Acts::MaterialSlab& slab,
                                         const RelativisticQuantities& rq) {
  // return early in case of vacuum or zero thickness
  if (not slab) {
    return 0.0f;
  }

  const float Ne = slab.material().molarElectronDensity();
  const float thickness = slab.thickness();
  // the Landau-Vavilov fwhm is 4*eps (see RPP2018 fig. 33.7)
  return 4 * computeEpsilon(Ne, thickness, rq);
}

}  // namespace detail

}  // namespace

float Acts::computeEnergyLossBethe(const MaterialSlab& slab, float m,
                                   float qOverP, float absQ) {
  // return early in case of vacuum or zero thickness
  if (not slab) {
    return 0.0f;
  }

  const RelativisticQuantities rq{m, qOverP, absQ};
  const float I = slab.material().meanExcitationEnergy();
  const float Ne = slab.material().molarElectronDensity();
  const float thickness = slab.thickness();
  const float eps = computeEpsilon(Ne, thickness, rq);
  const float dhalf = computeDeltaHalf(I, Ne, rq);
  const float u = computeMassTerm(Me, rq);
  const float wmax = computeWMax(m, rq);
  // uses RPP2018 eq. 33.5 scaled from mass stopping power to linear stopping
  // power and multiplied with the material thickness to get a total energy loss
  // instead of an energy loss per length.
  // the required modification only change the prefactor which becomes
  // identical to the prefactor epsilon for the most probable value.
  const float running =
      std::log(u / I) + std::log(wmax / I) - 2.0f * rq.beta2 - 2.0f * dhalf;
  return eps * running;
}

float Acts::deriveEnergyLossBetheQOverP(const MaterialSlab& slab, float m,
                                        float qOverP, float absQ) {
  // return early in case of vacuum or zero thickness
  if (not slab) {
    return 0.0f;
  }

  const RelativisticQuantities rq{m, qOverP, absQ};
  const float I = slab.material().meanExcitationEnergy();
  const float Ne = slab.material().molarElectronDensity();
  const float thickness = slab.thickness();
  const float eps = computeEpsilon(Ne, thickness, rq);
  const float dhalf = computeDeltaHalf(I, Ne, rq);
  const float u = computeMassTerm(Me, rq);
  const float wmax = computeWMax(m, rq);
  // original equation is of the form
  //
  //     eps * (log(u/I) + log(wmax/I) - 2 * beta² - delta)
  //     = eps * (log(u) + log(wmax) - 2 * log(I) - 2 * beta² - delta)
  //
  // with the resulting derivative as
  //
  //     d(eps) * (log(u/I) + log(wmax/I) - 2 * beta² - delta)
  //     + eps * (d(u)/(u) + d(wmax)/(wmax) - 2 * d(beta²) - d(delta))
  //
  // where we can use d(eps) = eps * (d(eps)/eps) for further simplification.
  const float logDerEps = logDeriveEpsilon(qOverP, rq);
  const float derDHalf = deriveDeltaHalf(qOverP, rq);
  const float logDerU = logDeriveMassTerm(qOverP);
  const float logDerWmax = logDeriveWMax(m, qOverP, rq);
  const float derBeta2 = deriveBeta2(qOverP, rq);
  const float rel = logDerEps * (std::log(u / I) + std::log(wmax / I) -
                                 2.0f * rq.beta2 - 2.0f * dhalf) +
                    logDerU + logDerWmax - 2.0f * derBeta2 - 2.0f * derDHalf;
  return eps * rel;
}

float Acts::computeEnergyLossLandau(const MaterialSlab& slab, float m,
                                    float qOverP, float absQ) {
  // return early in case of vacuum or zero thickness
  if (not slab) {
    return 0.0f;
  }

  const RelativisticQuantities rq{m, qOverP, absQ};
  const float I = slab.material().meanExcitationEnergy();
  const float Ne = slab.material().molarElectronDensity();
  const float thickness = slab.thickness();
  const float eps = computeEpsilon(Ne, thickness, rq);
  const float dhalf = computeDeltaHalf(I, Ne, rq);
  const float u = computeMassTerm(Me, rq);
  // uses RPP2018 eq. 33.12
  const float running =
      std::log(u / I) + std::log(eps / I) + 0.2f - rq.beta2 - 2 * dhalf;
  return eps * running;
}

float Acts::deriveEnergyLossLandauQOverP(const MaterialSlab& slab, float m,
                                         float qOverP, float absQ) {
  // return early in case of vacuum or zero thickness
  if (not slab) {
    return 0.0f;
  }

  const RelativisticQuantities rq{m, qOverP, absQ};
  const float I = slab.material().meanExcitationEnergy();
  const float Ne = slab.material().molarElectronDensity();
  const float thickness = slab.thickness();
  const float eps = computeEpsilon(Ne, thickness, rq);
  const float dhalf = computeDeltaHalf(I, Ne, rq);
  const float t = computeMassTerm(Me, rq);
  // original equation is of the form
  //
  //     eps * (log(t/I) - log(eps/I) - 0.2 - beta² - delta)
  //     = eps * (log(t) - log(eps) - 2*log(I) - 0.2 - beta² - delta)
  //
  // with the resulting derivative as
  //
  //     d(eps) * (log(t/I) - log(eps/I) - 0.2 - beta² - delta)
  //     + eps * (d(t)/t + d(eps)/eps - d(beta²) - 2*d(delta/2))
  //
  // where we can use d(eps) = eps * (d(eps)/eps) for further simplification.
  const float logDerEps = logDeriveEpsilon(qOverP, rq);
  const float derDHalf = deriveDeltaHalf(qOverP, rq);
  const float logDerT = logDeriveMassTerm(qOverP);
  const float derBeta2 = deriveBeta2(qOverP, rq);
  const float rel = logDerEps * (std::log(t / I) + std::log(eps / I) - 0.2f -
                                 rq.beta2 - 2 * dhalf) +
                    logDerT + logDerEps - derBeta2 - 2 * derDHalf;
  return eps * rel;
}

float Acts::computeEnergyLossLandauSigma(const MaterialSlab& slab, float m,
                                         float qOverP, float absQ) {
  // return early in case of vacuum or zero thickness
  if (not slab) {
    return 0.0f;
  }

  const RelativisticQuantities rq{m, qOverP, absQ};
  const float Ne = slab.material().molarElectronDensity();
  const float thickness = slab.thickness();
  // the Landau-Vavilov fwhm is 4*eps (see RPP2018 fig. 33.7)
  const float fwhm = 4 * computeEpsilon(Ne, thickness, rq);
  return convertLandauFwhmToGaussianSigma(fwhm);
}

float Acts::computeEnergyLossLandauFwhm(const MaterialSlab& slab, float m,
                                        float qOverP, float absQ) {
  const RelativisticQuantities rq{m, qOverP, absQ};
  return detail::computeEnergyLossLandauFwhm(slab, rq);
}

float Acts::computeEnergyLossLandauSigmaQOverP(const MaterialSlab& slab,
                                               float m, float qOverP,
                                               float absQ) {
  const RelativisticQuantities rq{m, qOverP, absQ};
  const float fwhm = detail::computeEnergyLossLandauFwhm(slab, rq);
  const float sigmaE = convertLandauFwhmToGaussianSigma(fwhm);
  //  var(q/p) = (d(q/p)/dE)² * var(E)
  // d(q/p)/dE = d/dE (q/sqrt(E²-m²))
  //           = q * -(1/2) * 1/p³ * 2E
  //           = -q/p² E/p = -(q/p)² * 1/(q*beta) = -(q/p)² * (q/beta) / q²
  //  var(q/p) = (q/p)^4 * (q/beta)² * (1/q)^4 * var(E)
  //           = (1/p)^4 * (q/beta)² * var(E)
  // do not need to care about the sign since it is only used squared
  const float pInv = qOverP / absQ;
  const float qOverBeta = std::sqrt(rq.q2OverBeta2);
  return qOverBeta * pInv * pInv * sigmaE;
}

namespace {

inline float computeBremsstrahlungLossMean(float mass, float energy) {
  return energy * (Me / mass) * (Me / mass);
}

inline float deriveBremsstrahlungLossMeanE(float mass) {
  return (Me / mass) * (Me / mass);
}

constexpr float MuonHighLowThreshold = 1_TeV;
// [low0 / X0] = MeV / mm -> [low0] = MeV
constexpr double MuonLow0 = -0.5345_MeV;
// [low1 * E / X0] = MeV / mm -> [low1] = 1
constexpr double MuonLow1 = 6.803e-5;
// [low2 * E^2 / X0] = MeV / mm -> [low2] = 1/MeV
constexpr double MuonLow2 = 2.278e-11 / 1_MeV;
// [low3 * E^3 / X0] = MeV / mm -> [low3] = 1/MeV^2
constexpr double MuonLow3 = -9.899e-18 / (1_MeV * 1_MeV);
// units are the same as low0
constexpr double MuonHigh0 = -2.986_MeV;
// units are the same as low1
constexpr double MuonHigh1 = 9.253e-5;

inline float computeMuonDirectPairPhotoNuclearLossMean(double energy) {
  if (energy < MuonHighLowThreshold) {
    return MuonLow0 +
           (MuonLow1 + (MuonLow2 + MuonLow3 * energy) * energy) * energy;
  } else {
    return MuonHigh0 + MuonHigh1 * energy;
  }
}

inline float deriveMuonDirectPairPhotoNuclearLossMeanE(double energy) {
  if (energy < MuonHighLowThreshold) {
    return MuonLow1 + (2 * MuonLow2 + 3 * MuonLow3 * energy) * energy;
  } else {
    return MuonHigh1;
  }
}

}  // namespace

float Acts::computeEnergyLossRadiative(const MaterialSlab& slab,
                                       PdgParticle absPdg, float m,
                                       float qOverP, float absQ) {
  assert((absPdg == Acts::makeAbsolutePdgParticle(absPdg)) &&
         "pdg is not absolute");

  // return early in case of vacuum or zero thickness
  if (not slab) {
    return 0.0f;
  }

  // relative radiation length
  const float x = slab.thicknessInX0();
  // particle momentum and energy
  // do not need to care about the sign since it is only used squared
  const float momentum = absQ / qOverP;
  const float energy = std::hypot(m, momentum);

  float dEdx = computeBremsstrahlungLossMean(m, energy);

  // muon- or muon+
  // TODO magic number 8_GeV
  if ((absPdg == PdgParticle::eMuon) and (8_GeV < energy)) {
    dEdx += computeMuonDirectPairPhotoNuclearLossMean(energy);
  }
  // scale from energy loss per unit radiation length to total energy
  return dEdx * x;
}

float Acts::deriveEnergyLossRadiativeQOverP(const MaterialSlab& slab,
                                            PdgParticle absPdg, float m,
                                            float qOverP, float absQ) {
  assert((absPdg == Acts::makeAbsolutePdgParticle(absPdg)) &&
         "pdg is not absolute");

  // return early in case of vacuum or zero thickness
  if (not slab) {
    return 0.0f;
  }

  // relative radiation length
  const float x = slab.thicknessInX0();
  // particle momentum and energy
  // do not need to care about the sign since it is only used squared
  const float momentum = absQ / qOverP;
  const float energy = std::hypot(m, momentum);

  // compute derivative w/ respect to energy.
  float derE = deriveBremsstrahlungLossMeanE(m);

  // muon- or muon+
  // TODO magic number 8_GeV
  if ((absPdg == PdgParticle::eMuon) and (8_GeV < energy)) {
    derE += deriveMuonDirectPairPhotoNuclearLossMeanE(energy);
  }
  // compute derivative w/ respect to q/p by using the chain rule
  //     df(E)/d(q/p) = df(E)/dE dE/d(q/p)
  // with
  //     E = sqrt(m² + p²) = sqrt(m² + q²/(q/p)²)
  // and the resulting derivative
  //     dE/d(q/p) = -q² / ((q/p)³ * E)
  const float derQOverP = -(absQ * absQ) / (qOverP * qOverP * qOverP * energy);
  return derE * derQOverP * x;
}

float Acts::computeEnergyLossMean(const MaterialSlab& slab, PdgParticle absPdg,
                                  float m, float qOverP, float absQ) {
  return computeEnergyLossBethe(slab, m, qOverP, absQ) +
         computeEnergyLossRadiative(slab, absPdg, m, qOverP, absQ);
}

float Acts::deriveEnergyLossMeanQOverP(const MaterialSlab& slab,
                                       PdgParticle absPdg, float m,
                                       float qOverP, float absQ) {
  return deriveEnergyLossBetheQOverP(slab, m, qOverP, absQ) +
         deriveEnergyLossRadiativeQOverP(slab, absPdg, m, qOverP, absQ);
}

float Acts::computeEnergyLossMode(const MaterialSlab& slab, PdgParticle absPdg,
                                  float m, float qOverP, float absQ) {
  // see ATL-SOFT-PUB-2008-003 section 3 for the relative fractions
  // TODO this is inconsistent with the text of the note
  return 0.9f * computeEnergyLossLandau(slab, m, qOverP, absQ) +
         0.15f * computeEnergyLossRadiative(slab, absPdg, m, qOverP, absQ);
}

float Acts::deriveEnergyLossModeQOverP(const MaterialSlab& slab,
                                       PdgParticle absPdg, float m,
                                       float qOverP, float absQ) {
  // see ATL-SOFT-PUB-2008-003 section 3 for the relative fractions
  // TODO this is inconsistent with the text of the note
  return 0.9f * deriveEnergyLossLandauQOverP(slab, m, qOverP, absQ) +
         0.15f * deriveEnergyLossRadiativeQOverP(slab, absPdg, m, qOverP, absQ);
}

namespace {

inline float theta0Highland(float xOverX0, float momentumInv,
                            float q2OverBeta2) {
  // RPP2018 eq. 33.15 (treats beta and q² consistently)
  const float t = std::sqrt(xOverX0 * q2OverBeta2);
  // log((x/X0) * (q²/beta²)) = log((sqrt(x/X0) * (q/beta))²)
  //                          = 2 * log(sqrt(x/X0) * (q/beta))
  return 13.6_MeV * momentumInv * t * (1.0f + 0.038f * 2 * std::log(t));
}

inline float theta0RossiGreisen(float xOverX0, float momentumInv,
                                float q2OverBeta2) {
  // TODO add source paper/ resource
  const float t = std::sqrt(xOverX0 * q2OverBeta2);
  return 17.5_MeV * momentumInv * t *
         (1.0f + 0.125f * std::log10(10.0f * xOverX0));
}

}  // namespace

float Acts::computeMultipleScatteringTheta0(const MaterialSlab& slab,
                                            PdgParticle absPdg, float m,
                                            float qOverP, float absQ) {
  assert((absPdg == Acts::makeAbsolutePdgParticle(absPdg)) &&
         "pdg is not absolute");

  // return early in case of vacuum or zero thickness
  if (not slab) {
    return 0.0f;
  }

  // relative radiation length
  const float xOverX0 = slab.thicknessInX0();
  // 1/p = q/(pq) = (q/p)/q
  const float momentumInv = std::abs(qOverP / absQ);
  // q²/beta²; a smart compiler should be able to remove the unused computations
  const float q2OverBeta2 = RelativisticQuantities(m, qOverP, absQ).q2OverBeta2;

  // electron or positron
  if (absPdg == PdgParticle::eElectron) {
    return theta0RossiGreisen(xOverX0, momentumInv, q2OverBeta2);
  } else {
    return theta0Highland(xOverX0, momentumInv, q2OverBeta2);
  }
}