diff --git a/Kassiopeia/Bindings/CMakeLists.txt b/Kassiopeia/Bindings/CMakeLists.txt index daea9343f..2ab7c85d0 100644 --- a/Kassiopeia/Bindings/CMakeLists.txt +++ b/Kassiopeia/Bindings/CMakeLists.txt @@ -119,6 +119,7 @@ set( BINDINGS_HEADER_FILES Interactions/Include/KSIntCalculatorHydrogenBuilder.h Interactions/Include/KSIntCalculatorIonBuilder.h Interactions/Include/KSIntCalculatorArgonBuilder.h + Interactions/Include/KSIntCalculatorMottBuilder.h Interactions/Include/KSIntCalculatorKESSBuilder.h Interactions/Include/KSIntScatteringBuilder.h Interactions/Include/KSIntSpinFlipBuilder.h @@ -346,6 +347,7 @@ set( BINDINGS_SOURCE_FILES Interactions/Source/KSIntCalculatorHydrogenBuilder.cxx Interactions/Source/KSIntCalculatorIonBuilder.cxx Interactions/Source/KSIntCalculatorArgonBuilder.cxx + Interactions/Source/KSIntCalculatorMottBuilder.cxx Interactions/Source/KSIntCalculatorKESSBuilder.cxx Interactions/Source/KSIntScatteringBuilder.cxx Interactions/Source/KSIntSpinFlipBuilder.cxx diff --git a/Kassiopeia/Bindings/Interactions/Include/KSIntCalculatorMottBuilder.h b/Kassiopeia/Bindings/Interactions/Include/KSIntCalculatorMottBuilder.h new file mode 100644 index 000000000..e843f285c --- /dev/null +++ b/Kassiopeia/Bindings/Interactions/Include/KSIntCalculatorMottBuilder.h @@ -0,0 +1,62 @@ +#ifndef Kassiopeia_KSIntCalculatorMottBuilder_h_ +#define Kassiopeia_KSIntCalculatorMottBuilder_h_ + +#include "KComplexElement.hh" +#include "KSIntCalculatorMott.h" + +using namespace Kassiopeia; +namespace katrin +{ + +typedef KComplexElement KSIntCalculatorMottBuilder; + +template<> inline bool KSIntCalculatorMottBuilder::AddAttribute(KContainer* aContainer) +{ + if (aContainer->GetName() == "name") { + aContainer->CopyTo(fObject, &KNamed::SetName); + return true; + } + if (aContainer->GetName() == "theta_min") { + aContainer->CopyTo(fObject, &KSIntCalculatorMott::SetThetaMin); + + if ((fObject->GetThetaMin() <= 0.0)) { + initmsg(eError) << "\"" << fObject->GetThetaMin() + << R"(" is an invalid lower bound for Mott scattering! Change to a value > 0)" << eom; + + return false; + } + else { + return true; + } + } + if (aContainer->GetName() == "theta_max") { + aContainer->CopyTo(fObject, &KSIntCalculatorMott::SetThetaMax); + if ((fObject->GetThetaMax() > 180.0)) { + initmsg(eError) << "\"" << fObject->GetThetaMax() + << R"(" is an invalid upper bound for Mott scattering! Change to a value < 180)" << eom; + + return false; + } + else { + return true; + } + } + if (aContainer->GetName() == "nucleus") { + aContainer->CopyTo(fObject, &KSIntCalculatorMott::SetNucleus); + + if ((fObject->GetNucleus().compare("He-4") != 0) && (fObject->GetNucleus().compare("Ne-20") != 0)) { + initmsg(eError) << "\"" << fObject->GetNucleus() + << R"(" is not available for Mott scattering! Available gases: "He-4", "Ne-20")" << eom; + + return false; + } + else { + return true; + } + } + return false; +} + +} // namespace katrin + +#endif diff --git a/Kassiopeia/Bindings/Interactions/Source/KSIntCalculatorMottBuilder.cxx b/Kassiopeia/Bindings/Interactions/Source/KSIntCalculatorMottBuilder.cxx new file mode 100644 index 000000000..58891e168 --- /dev/null +++ b/Kassiopeia/Bindings/Interactions/Source/KSIntCalculatorMottBuilder.cxx @@ -0,0 +1,20 @@ +#include "KSIntCalculatorMottBuilder.h" + +#include "KSRootBuilder.h" + +using namespace Kassiopeia; +using namespace std; + +namespace katrin +{ + +template<> KSIntCalculatorMottBuilder::~KComplexElement() = default; + +STATICINT sKSIntCalculatorMottStructure = KSIntCalculatorMottBuilder::Attribute("name") + + KSIntCalculatorMottBuilder::Attribute("theta_min") + + KSIntCalculatorMottBuilder::Attribute("theta_max") + + KSIntCalculatorMottBuilder::Attribute("nucleus"); + +STATICINT sToolboxKSIntCalculatorMott = + KSRootBuilder::ComplexElement("ksint_calculator_mott"); +} // namespace katrin diff --git a/Kassiopeia/Interactions/CMakeLists.txt b/Kassiopeia/Interactions/CMakeLists.txt index d423ccc92..cb9ee5b9c 100644 --- a/Kassiopeia/Interactions/CMakeLists.txt +++ b/Kassiopeia/Interactions/CMakeLists.txt @@ -24,6 +24,7 @@ set( INTERACTIONS_HEADER_BASENAMES KSIntCalculatorIon.h KSIntCalculatorTritium.h KSIntCalculatorArgon.h + KSIntCalculatorMott.h KSIntSurfaceSpecular.h KSIntSurfaceUCN.h @@ -78,6 +79,7 @@ set( INTERACTIONS_SOURCE_BASENAMES KSIntCalculatorIon.cxx KSIntCalculatorTritium.cxx KSIntCalculatorArgon.cxx + KSIntCalculatorMott.cxx KSIntSurfaceSpecular.cxx KSIntSurfaceUCN.cxx diff --git a/Kassiopeia/Interactions/Include/KSIntCalculatorMott.h b/Kassiopeia/Interactions/Include/KSIntCalculatorMott.h new file mode 100644 index 000000000..3bdd9c389 --- /dev/null +++ b/Kassiopeia/Interactions/Include/KSIntCalculatorMott.h @@ -0,0 +1,127 @@ +#ifndef Kassiopeia_KSIntCalculatorMott_h_ +#define Kassiopeia_KSIntCalculatorMott_h_ + +#include "KField.h" +#include "KSIntCalculator.h" + +/* + * KSintCalculatorMott.h + * + * Date: August 22, 2022 + * Author: Junior Ivan Peña (juniorpe) + */ + +namespace Kassiopeia +{ + +/* + * The xml configuration for using this calculator is as follows: + * + * where [lower_bound] is the lower limit, in degrees, on the range of allowed scattering + * angles and [upper_bound] is the upper limit. A theta_min of 0 is not allowed since the + * Mott Cross Section has a singularity there. + */ +class KSIntCalculatorMott : public KSComponentTemplate +{ + public: + KSIntCalculatorMott(); + KSIntCalculatorMott(const KSIntCalculatorMott& aCopy); + KSIntCalculatorMott* Clone() const override; + ~KSIntCalculatorMott() override; + + public: + void CalculateCrossSection(const KSParticle& anInitialParticle, double& aCrossSection) override; + void ExecuteInteraction(const KSParticle& anInitialParticle, KSParticle& aFinalParticle, + KSParticleQueue& aSecondaries) override; + + public: + ; + K_SET_GET(double, ThetaMin); // radians + ; + K_SET_GET(double, ThetaMax); // radians + ; + K_SET_GET(std::string, Nucleus); + + public: + /** + * \brief Returns scattering angle sampled from Mott Differential Cross Seciton for a given + * incoming electron energy. Sampling is done using using inverse transform sampling of + * Rutherford Differential Cross Section, then rejection sampling with a rescaled Ruth. Diff. X-Sec + * + * \parameter electron's initial kinetic energy + * + * \return Scatter angle + */ + double GetTheta(const double& anEnergy); + + /** + * \brief Calculates energy loss after electron scatters using equation (1.51) in: + * C. Leroy and P.G. Rancoita, Principles of Radiation Interaction in Matter and Detection, + * 2nd Edition, World Scientific (Singapore) 2009. + * + * \parameters electron's initial kinetic energy, scattering angle + * + * \return Electron's energy loss + */ + double GetEnergyLoss(const double& anEnergy, const double& aTheta); + /** + * \parameter Electron's initial kinetic energy + * + * \return Electron's natural velocity beta = v/c + */ + double Beta(double const E0); + /** + * \brief Takes in the electron's initial energy and calculates the coefficients used in the + * RMott calculation for He corresponding to equation (25) in: + * M.J Boschini, C. Consolandi, M. Gervasi, et al., J. Radiat. Phys. Chem. 90 (2013) 36-66. + * + * \parameter Electron's initial kinetic energy + * + * \return Vector containing 5 calculated coeffecients of a + */ + std::vector RMott_coeffs(double const E0); + /** + * \brief Returns the Mott Differential Cross Section(from equation (3) and (24) in: M.J Boschini, + * C. Consolandi, M. Gervasi, et al., J. Radiat. Phys. Chem. 90 (2013) 36-66.) Used for rejection + * sampling in GetTheta function + * + * \parameter Electron's initial kinetic energy and scatter angle + * + * \return Mott Differential Cross Section + * + */ + double MDCS(double theta, const double E0); + /** + * \brief Calculates the analytical integral result of the Mott Differential Cross Section(from + * equation (3) and (24) in: M.J Boschini, C. Consolandi, M. Gervasi, et al., J. Radiat. Phys. Chem. 90 (2013) 36-66.) + * within the bounds [fThetaMin, theta] for a given electron initial kinetic energy and some + * integral upper bound theta + * + * \parameter Electron's initial kinetic energy, integral upper bound theta + * + * \return Integral of Mott Differential Cross Section from fThetaMin to theta + * + */ + double MTCS(double theta, const double E0); + /** + * \brief Used for rejection sampling in GetTheta function + * + * \parameter Scatter Angle + * + * \return Rutherford Differential Cross Section normalized in range [fThetaMin, fThetaMax] + * + */ + double Normalized_RDCS(double theta); + /** + * \brief Used for inverse transform sampling for obtaining scatter angle in GetTheta function + * + * \return Inverse of normalized Rutherford Differential Cross Section in range [fThetaMin, fThetaMax] + * + */ + double Normalized_RTCSInverse(double x); + +}; + +} // namespace Kassiopeia + +#endif diff --git a/Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx b/Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx new file mode 100644 index 000000000..89c38ff15 --- /dev/null +++ b/Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx @@ -0,0 +1,224 @@ +#include "KSIntCalculatorMott.h" +#include +#include + +#include "KRandom.h" +using katrin::KRandom; + +#include "KConst.h" +#include "KStringUtils.h" +using katrin::KStringUtils; + +#include "KThreeVector.hh" +using katrin::KThreeVector; + +namespace Kassiopeia +{ + +KSIntCalculatorMott::KSIntCalculatorMott() : fThetaMin(0.), fThetaMax(0.), fNucleus("") {} +KSIntCalculatorMott::KSIntCalculatorMott(const KSIntCalculatorMott& aCopy) : + KSComponent(aCopy), + fThetaMin(aCopy.fThetaMin), + fThetaMax(aCopy.fThetaMax), + fNucleus(aCopy.fNucleus) +{} +KSIntCalculatorMott* KSIntCalculatorMott::Clone() const +{ + return new KSIntCalculatorMott(*this); +} +KSIntCalculatorMott::~KSIntCalculatorMott() = default; + +void KSIntCalculatorMott::CalculateCrossSection(const KSParticle& anInitialParticle, double& aCrossSection) +{ + + double tInitialKineticEnergy = anInitialParticle.GetKineticEnergy_eV(); + + double fCrossSection = (MTCS(fThetaMax, tInitialKineticEnergy) - MTCS(fThetaMin, tInitialKineticEnergy)); + + aCrossSection = fCrossSection; + return; +} +void KSIntCalculatorMott::ExecuteInteraction(const KSParticle& anInitialParticle, KSParticle& aFinalParticle, + KSParticleQueue& /*aSecondaries*/) +{ + double tTime = anInitialParticle.GetTime(); + double tInitialKineticEnergy = anInitialParticle.GetKineticEnergy_eV(); + KThreeVector tPosition = anInitialParticle.GetPosition(); + KThreeVector tMomentum = anInitialParticle.GetMomentum(); + + double tTheta = GetTheta(tInitialKineticEnergy); + double tLostKineticEnergy = GetEnergyLoss(tInitialKineticEnergy, tTheta); + + + double tPhi = KRandom::GetInstance().Uniform(0., 2. * katrin::KConst::Pi()); + + KThreeVector tOrthogonalOne = tMomentum.Orthogonal(); + KThreeVector tOrthogonalTwo = tMomentum.Cross(tOrthogonalOne); + KThreeVector tFinalDirection = + tMomentum.Magnitude() * + (sin(tTheta* katrin::KConst::Pi() / 180) * (cos(tPhi) * tOrthogonalOne.Unit() + sin(tPhi) * tOrthogonalTwo.Unit()) + + cos(tTheta* katrin::KConst::Pi() / 180) * tMomentum.Unit()); + + + aFinalParticle.SetTime(tTime); + aFinalParticle.SetPosition(tPosition); + aFinalParticle.SetMomentum(tFinalDirection); + aFinalParticle.SetKineticEnergy_eV(tInitialKineticEnergy - tLostKineticEnergy); + aFinalParticle.SetLabel(GetName()); + + fStepNInteractions = 1; + fStepEnergyLoss = tLostKineticEnergy; + fStepAngularChange = tTheta; + + return; +} + +double KSIntCalculatorMott::GetTheta(const double& anEnergy){ + + double tTheta; + double resolution; + double resolution_factor = 10; // arbitrarily chosen for o(1ms) computation time for He; computation time increases linearly with resolution factor; resolution factors under 1 yields incorrect results + + + resolution = ceil((fThetaMax - fThetaMin) * resolution_factor); + + + // Initializing array for possible theta values + std::vector theta_arr; + for (int i = 0; i <= int(resolution); ++i) { + theta_arr.push_back(fThetaMin + i * (fThetaMax - fThetaMin) / resolution); + } + + // Calculating rescaling factor for rejection sampling + std::vector k_values; + for (double theta : theta_arr) { + k_values.push_back(MDCS(theta, anEnergy) / Normalized_RDCS(theta)); + } + + // Find the maximum value of k + double k = *std::max_element(k_values.begin(), k_values.end()); + + + + // Rejection sampling taking initial sample from Rutherford Diff. X-Sec + while (true) { + double sample = Normalized_RTCSInverse(KRandom::GetInstance().Uniform(0.0, 1.0, false, true)); + double uniform_sample = KRandom::GetInstance().Uniform(0.0, 1.0, false, true); + if (MDCS(sample, anEnergy) / (k * Normalized_RDCS(sample)) > uniform_sample) { + tTheta = sample; + break; + } + } + + return tTheta; + +} + +double KSIntCalculatorMott::GetEnergyLoss(const double& anEnergy, const double& aTheta){ + double M, me, p, anELoss; + + double amu = 0.0; + + if (fNucleus == "He-4") { + amu = 4.0026; // mass in amu from W.J. Huang et al 2021 Chinese Phys. C 45 030002 + } else if (fNucleus == "Ne-20") { + amu = 19.9924; // mass in amu from W.J. Huang et al 2021 Chinese Phys. C 45 030002 + } + + M = amu * katrin::KConst::AtomicMassUnit_eV(); // eV/c^2 + me = katrin::KConst::M_el_eV(); // electron mass eV/c^2 + + p = sqrt(anEnergy * (anEnergy + 2 * me * pow(katrin::KConst::C(), 2))) / katrin::KConst::C(); + + anELoss = 2 * pow(p, 2) * M / (pow(me, 2) + pow(M, 2) + 2 * M * sqrt( pow((p/katrin::KConst::C()), 2) + pow(me, 2))) * (1 - cos(aTheta * katrin::KConst::Pi() / 180)); + + return anELoss; +} + + +double KSIntCalculatorMott::Beta(double const E0) { + return sqrt( pow(E0, 2) + 2 * E0 * katrin::KConst::M_el_eV()) / (E0 + katrin::KConst::M_el_eV()); +} + +std::vector KSIntCalculatorMott::RMott_coeffs(double const E0) { + + std::vector a(6, 0.0); // Initialize a with 6 zeros, the last entry is not related to coefficient calculation it is Z for the nucleus of choice + std::vector> b(5, std::vector(6)); + + if (KStringUtils::IStartsWith(fNucleus, "He-")) { + a[5] = 2; // Charge Z + b = {{ 1.0 , 3.76476e-8, -3.05313e-7, -3.27422e-7, 2.44235e-6, 4.08754e-6}, + { 2.35767e-2, 3.24642e-2, -6.37269e-4, -7.69160e-4, 5.28004e-3, 9.45642e-3}, + {-2.73743e-1, -7.40767e-1, -4.98195e-1, 1.74337e-3, -1.25798e-2, -2.24046e-2}, + {-7.79128e-4, -4.14495e-4, -1.62657e-3, -1.37286e-3, 1.04319e-2, 1.83488e-2}, + { 2.02855e-4, 1.94598e-6, 4.30102e-4, 4.32180e-4, -3.31526e-3, -5.81788e-3}}; + } else if (KStringUtils::IStartsWith(fNucleus, "Ne-")) { + a[5] = 10; // Charge Z + b = {{ 9.99997e-1, -1.87404e-7, 3.10276e-5, 5.20000e-5, 2.98132e-4, -5.19259e-4}, + { 1.20783e-1, 1.66407e-1, 1.06608e-2, 6.48772e-3, -1.53031e-3, -7.59354e-2}, + {-3.15222e-1, -8.28793e-1, -6.05740e-1, -1.47812e-1, 1.15760 , 1.58565 }, + {-2.89055e-2, -9.08096e-4, 1.21467e-1, -1.77575e-1, -1.29110 , -1.90333 }, + { 7.65342e-3, -8.85417e-4, -3.89092e-2, -5.44040e-2, 3.93087e-1, 5.79439e-1}}; + } + + for (int j = 0; j < 5; ++j) { + a[j] = 0.0; + for (int k = 0; k < 6; ++k) { + a[j] += b[j][k] * pow(Beta(E0) - 0.7181287, k); + } + } + + return a; +} + +double KSIntCalculatorMott::MDCS(double theta, const double E0) { + + double r_e = 2.8179403227 * pow(10, -15); // classical electron radius + + std::vector a = RMott_coeffs(E0); + + double Z = a[5]; + + return pow(Z, 2) * pow(r_e, 2) * ( (1 - pow(Beta(E0), 2))/(pow(Beta(E0), 4)) ) * ((a[0] * pow(1 - cos(theta * katrin::KConst::Pi() / 180), -2)) + (a[1] * pow(sqrt(1 - cos(theta * katrin::KConst::Pi() / 180)), -3)) + (a[2] * pow((1 - cos(theta * katrin::KConst::Pi() / 180)), (-1))) + (a[3] * pow(sqrt(1 - cos(theta * katrin::KConst::Pi() / 180)), -1)) + a[4] ); +} + +double KSIntCalculatorMott::MTCS(double theta, const double E0) { + + double r_e = 2.8179403227 * pow(10, -15); // classical electron radius + + std::vector a = RMott_coeffs(E0); + + double Z = a[5]; + + + return 2 * katrin::KConst::Pi() * pow(Z, 2) * pow(r_e, 2) * + ((1 - pow(Beta(E0), 2)) / pow(Beta(E0), 4)) * + ((a[0] * pow((1 - cos(theta * katrin::KConst::Pi() / 180)), -1) / (-1)) + + (a[1] * pow((1 - cos(theta * katrin::KConst::Pi() / 180)), (1.0 / 2) - 1) / ((1.0 / 2) - 1)) + + (a[2] * log(1 - cos(theta * katrin::KConst::Pi() / 180))) + + (a[3] * pow((1 - cos(theta * katrin::KConst::Pi() / 180)), (3.0 / 2) - 1) / ((3.0 / 2) - 1)) + + (a[4] * -1 * cos(theta * katrin::KConst::Pi() / 180))) + - 2 * katrin::KConst::Pi() * pow(Z, 2) * pow(r_e, 2) * + ((1 - pow(Beta(E0), 2)) / pow(Beta(E0), 4)) * + ((a[0] * pow((1 - cos(fThetaMin * katrin::KConst::Pi() / 180)), -1) / (-1)) + + (a[1] * pow((1 - cos(fThetaMin * katrin::KConst::Pi() / 180)), (1.0 / 2) - 1) / ((1.0 / 2) - 1)) + + (a[2] * log(1 - cos(fThetaMin * katrin::KConst::Pi() / 180))) + + (a[3] * pow((1 - cos(fThetaMin * katrin::KConst::Pi() / 180)), (3.0 / 2) - 1) / ((3.0 / 2) - 1)) + + (a[4] * -1 * cos(fThetaMin * katrin::KConst::Pi() / 180))); +} + +double KSIntCalculatorMott::Normalized_RDCS(double theta) { + + double C = 1 / ((1 / (1 - cos(fThetaMin* katrin::KConst::Pi() / 180))) - (1 / (1 - cos(fThetaMax* katrin::KConst::Pi() / 180)))); + + return C / pow((1 - cos(theta* katrin::KConst::Pi() / 180)), 2); +} + +double KSIntCalculatorMott::Normalized_RTCSInverse(double x) { + + double C = 1 / ((1 / (1 - cos(fThetaMin * katrin::KConst::Pi() / 180))) - (1 / (1 - cos(fThetaMax * katrin::KConst::Pi() / 180)))); + + return acos( 1 - (1 / ((1 / (1 - cos(fThetaMin * katrin::KConst::Pi() / 180))) - x / C)) ) * 180 / katrin::KConst::Pi(); +} + +} // namespace Kassiopeia diff --git a/Kassiopeia/XML/Complete/Interactions.xml b/Kassiopeia/XML/Complete/Interactions.xml index 0e92da518..5cd7cd4aa 100644 --- a/Kassiopeia/XML/Complete/Interactions.xml +++ b/Kassiopeia/XML/Complete/Interactions.xml @@ -112,4 +112,33 @@ only H_2 is available at present and is used by default. --> + + + + + + +