From bfae4037e68bb667b323d200bb2a7afee4d5e338 Mon Sep 17 00:00:00 2001 From: juniorpena <72262395+juniorpena@users.noreply.github.com> Date: Mon, 22 Aug 2022 16:55:44 -0400 Subject: [PATCH 01/22] Add files via upload --- .../Include/KSIntCalculatorMott.h | 87 +++++++++++++++++++ 1 file changed, 87 insertions(+) create mode 100644 Kassiopeia/Interactions/Include/KSIntCalculatorMott.h diff --git a/Kassiopeia/Interactions/Include/KSIntCalculatorMott.h b/Kassiopeia/Interactions/Include/KSIntCalculatorMott.h new file mode 100644 index 000000000..b92c59631 --- /dev/null +++ b/Kassiopeia/Interactions/Include/KSIntCalculatorMott.h @@ -0,0 +1,87 @@ +#ifndef Kassiopeia_KSIntCalculatorMott_h_ +#define Kassiopeia_KSIntCalculatorMott_h_ + +#include "KField.h" +#include "KSIntCalculator.h" +#include "TF1.h" + +/* + * KSintCalculatorMott.h + * + * Date: August 22, 2022 + * Author: Junior Peña (juniorpe) + */ + +namespace Kassiopeia +{ + +/* + * The xml configuration for using this calculator is as follows: + * + * where [lower_bound] is the lower limit, in radians, 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 + + public: + /** + * \brief Returns scattering angle sampled from Mott Differential Cross Seciton for a given + * incoming electron energy. Sampling is done using the GetRandom method for ROOT's TF1 class. + * + * \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); + /** + * \brief Initializes Mott Differential Cross Section given in: + * M.J Boschini, C. Consolandi, M. Gervasi, et al., J. Radiat. Phys. Chem. 90 (2013) 36-66. + * Also initializes Total Cross Section, where integration was done separately in mathematica + * and analytical form written in terms of constants given in publication above. The cross sections + * are initialized as ROOT TF1 objects. + * + */ + void InitializeMDCS(const double E0); + /** + * \brief Deinitializes Mott Differential Cross Section and Mott Total Cross Section + * + */ + void DeinitializeMDCS(); + + protected: + TF1* fMDCS; + TF1* fMTCS; +}; + +} // namespace Kassiopeia + +#endif From a1b0c2f87422301070cc507460eaa867d5f589e0 Mon Sep 17 00:00:00 2001 From: juniorpena <72262395+juniorpena@users.noreply.github.com> Date: Mon, 22 Aug 2022 16:57:22 -0400 Subject: [PATCH 02/22] Add files via upload --- .../Source/KSIntCalculatorMott.cxx | 175 ++++++++++++++++++ 1 file changed, 175 insertions(+) create mode 100644 Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx diff --git a/Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx b/Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx new file mode 100644 index 000000000..c9d2b2553 --- /dev/null +++ b/Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx @@ -0,0 +1,175 @@ +#include "KSIntCalculatorMott.h" + +#include "KRandom.h" +using katrin::KRandom; + +#include "KConst.h" + +using KGeoBag::KThreeVector; + +namespace Kassiopeia +{ + +KSIntCalculatorMott::KSIntCalculatorMott() : fThetaMin(0.), fThetaMax(0.), fMDCS(nullptr), fMTCS(nullptr) {} +KSIntCalculatorMott::KSIntCalculatorMott(const KSIntCalculatorMott& aCopy) : + KSComponent(aCopy), + fThetaMin(aCopy.fThetaMin), + fThetaMax(aCopy.fThetaMax), + fMDCS(nullptr), + fMTCS(nullptr) +{} +KSIntCalculatorMott* KSIntCalculatorMott::Clone() const +{ + return new KSIntCalculatorMott(*this); +} +KSIntCalculatorMott::~KSIntCalculatorMott() = default; + +void KSIntCalculatorMott::CalculateCrossSection(const KSParticle& anInitialParticle, double& aCrossSection) +{ + + double tInitialKineticEnergy = anInitialParticle.GetKineticEnergy_eV(); + InitializeMDCS(tInitialKineticEnergy); + + double fCrossSection = (fMTCS->Eval(fThetaMax) - fMTCS->Eval(fThetaMin)) * pow(10, -30); + + DeinitializeMDCS(); + + 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()); + tMomentum.SetPolarAngle(tTheta); + tMomentum.SetAzimuthalAngle(tPhi); + + aFinalParticle.SetTime(tTime); + aFinalParticle.SetPosition(tPosition); + aFinalParticle.SetMomentum(tMomentum); + aFinalParticle.SetKineticEnergy_eV(tInitialKineticEnergy - tLostKineticEnergy); + aFinalParticle.SetLabel(GetName()); + + fStepAngularChange = tTheta * 180. / katrin::KConst::Pi(); + + return; +} + +double KSIntCalculatorMott::GetTheta(const double& anEnergy){ + + double tTheta; + + InitializeMDCS(anEnergy); + + tTheta = fMDCS->GetRandom(fThetaMin, fThetaMax); + + DeinitializeMDCS(); + + return tTheta; + +} + +double KSIntCalculatorMott::GetEnergyLoss(const double& anEnergy, const double& aTheta){ + double M, me, p, anELoss; + + M = 4.0026 * 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)); + + return anELoss; +} + +void KSIntCalculatorMott::InitializeMDCS(const double E0) +{ + + double beta = TMath::Sqrt( pow(E0, 2) + 2 * E0 * katrin::KConst::M_el_eV()) / (E0 + katrin::KConst::M_el_eV()); + + /* Constants given in publication in header file. These are the constants for scattering off of Helium */ + static double b[5][6] = {{ 1. , 3.76476 * pow(10, -8), -3.05313 * pow(10, -7), -3.27422 * pow(10, -7), 2.44235 * pow(10, -6), 4.08754 * pow(10, -6)}, + { 2.35767 * pow(10, -2), 3.24642 * pow(10, -2), -6.37269 * pow(10, -4), -7.69160 * pow(10, -4), 5.28004 * pow(10, -3), 9.45642 * pow(10, -3)}, + { -2.73743 * pow(10, -1), -7.40767 * pow(10, -1), -4.98195 * pow(10, -1), 1.74337 * pow(10, -3), -1.25798 * pow(10, -2), -2.24046 * pow(10, -2)}, + { -7.79128 * pow(10, -4), -4.14495 * pow(10, -4), -1.62657 * pow(10, -3), -1.37286 * pow(10, -3), 1.04319 * pow(10, -2), 1.83488 * pow(10, -2)}, + { 2.02855 * pow(10, -4), 1.94598 * pow(10, -6), 4.30102 * pow(10, -4), 4.32180 * pow(10, -4), -3.31526 * pow(10, -3), -5.81788 * pow(10, -3)}}; + + double a[6] = {0, 0, 0, 0, 0, 0}; + + for (int j = 0; j < 5; j++){ + for (int k = 0; k < 6; k++){ + a[j] += b[j][k] * pow((beta - 0.7181287), k); + } + } + + /* ROOT TF1 takes in analytical formulas as strings. This is the conversion into strings */ + std::string Rmott; + std::ostringstream RmottStream; + + RmottStream << std::fixed; + RmottStream << std::setprecision(12); + + RmottStream << a[0] << " + " << a[1] << " * pow( ( 1 - cos(x) ) , 0.5) + " << a[2] << " * ( 1 - cos(x) ) + " << a[3] << " * pow( ( 1 - cos(x) ) , 1.5) + " << a[4] << " * pow( ( 1 - cos(x) ) , 2)"; + Rmott = RmottStream.str(); + + + + std::string RutherfordDCS; + std::ostringstream RutherfordDCSStream; + + RutherfordDCSStream << std::fixed; + RutherfordDCSStream << std::setprecision(12); + + RutherfordDCSStream << "pow((2 * 2.8179403227 ), 2) * ( (1 - pow(" << beta << ", 2) ) / ( pow(" << beta << ", 4) )) * (1 / pow( (1 - cos(x)), 2))"; + RutherfordDCS = RutherfordDCSStream.str(); + + + + std::string MottDCS = "(" + RutherfordDCS + ") * (" + Rmott + ")"; + + + std::string MottTCS; + std::ostringstream MottTCSStream; + std::ostringstream MottTCSStreamThetaIndepPrefactor; + + MottTCSStream << std::fixed; + MottTCSStream << std::setprecision(12); + + MottTCSStreamThetaIndepPrefactor << std::fixed; + MottTCSStreamThetaIndepPrefactor << std::setprecision(12); + + MottTCSStreamThetaIndepPrefactor << "pow((2 * 2.8179403227 ), 2) * ( (1 - pow(" << beta << ", 2) ) / ( pow(" << beta << ", 4) )) * "; + + MottTCSStream << "(2 * " << katrin::KConst::Pi() << " * (-4 * sqrt(2) * ( (" << a[1] << ") + (-1) * (" << a[2] + << ") * sqrt(( 1 - cos(x) )) * log( sin(x / 2) )) * pow(sin(x / 2), 4) + 4 * sqrt(2) * (2 * (" + << a[3] << ") + (" << a[4] << ") * sqrt( (1 - cos(x))) ) * pow(sin(x / 2), 6) - 2 * (" << a[0] + << ") * pow(sin( x / 2 ), 3) ) ) / ( pow((-1 + cos(x)), 2) * sin( x / 2 ))"; + + MottTCS = MottTCSStreamThetaIndepPrefactor.str() + MottTCSStream.str(); + + fMDCS = new TF1("function", MottDCS.c_str(), fThetaMin, fThetaMax); + fMTCS = new TF1("function", MottTCS.c_str(), fThetaMin, fThetaMax); + + return; +} + +void KSIntCalculatorMott::DeinitializeMDCS() +{ + delete fMDCS; + fMDCS = nullptr; + delete fMTCS; + fMTCS = nullptr; + return; +} + +} // namespace Kassiopeia From 83ce264ef3a16c1713eb5642ec1b763bda4a2a41 Mon Sep 17 00:00:00 2001 From: juniorpena <72262395+juniorpena@users.noreply.github.com> Date: Mon, 22 Aug 2022 16:58:56 -0400 Subject: [PATCH 03/22] Add files via upload --- .../Include/KSIntCalculatorMottBuilder.h | 32 +++++++++++++++++++ 1 file changed, 32 insertions(+) create mode 100644 Kassiopeia/Bindings/Interactions/Include/KSIntCalculatorMottBuilder.h diff --git a/Kassiopeia/Bindings/Interactions/Include/KSIntCalculatorMottBuilder.h b/Kassiopeia/Bindings/Interactions/Include/KSIntCalculatorMottBuilder.h new file mode 100644 index 000000000..6f6880f2c --- /dev/null +++ b/Kassiopeia/Bindings/Interactions/Include/KSIntCalculatorMottBuilder.h @@ -0,0 +1,32 @@ +#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); + return true; + } + if (aContainer->GetName() == "theta_max") { + aContainer->CopyTo(fObject, &KSIntCalculatorMott::SetThetaMax); + return true; + } + return false; +} + +} // namespace katrin + +#endif From 96f00f5f1b966a2daaf0850f08a2575fadf1754f Mon Sep 17 00:00:00 2001 From: juniorpena <72262395+juniorpena@users.noreply.github.com> Date: Mon, 22 Aug 2022 16:59:41 -0400 Subject: [PATCH 04/22] Add files via upload --- .../Source/KSIntCalculatorMottBuilder.cxx | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) create mode 100644 Kassiopeia/Bindings/Interactions/Source/KSIntCalculatorMottBuilder.cxx diff --git a/Kassiopeia/Bindings/Interactions/Source/KSIntCalculatorMottBuilder.cxx b/Kassiopeia/Bindings/Interactions/Source/KSIntCalculatorMottBuilder.cxx new file mode 100644 index 000000000..7590801d0 --- /dev/null +++ b/Kassiopeia/Bindings/Interactions/Source/KSIntCalculatorMottBuilder.cxx @@ -0,0 +1,19 @@ +#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"); + +STATICINT sToolboxKSIntCalculatorMott = + KSRootBuilder::ComplexElement("ksint_calculator_mott"); +} // namespace katrin From 70eadb548f7a0f2c987b375abd79e860215ad430 Mon Sep 17 00:00:00 2001 From: juniorpena <72262395+juniorpena@users.noreply.github.com> Date: Mon, 22 Aug 2022 17:01:24 -0400 Subject: [PATCH 05/22] Update CMakeLists.txt --- Kassiopeia/Interactions/CMakeLists.txt | 2 ++ 1 file changed, 2 insertions(+) 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 From 86820c78bedbec3aca760d84d60ac13a4afec36e Mon Sep 17 00:00:00 2001 From: juniorpena <72262395+juniorpena@users.noreply.github.com> Date: Mon, 22 Aug 2022 17:03:07 -0400 Subject: [PATCH 06/22] Update CMakeLists.txt --- Kassiopeia/Bindings/CMakeLists.txt | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Kassiopeia/Bindings/CMakeLists.txt b/Kassiopeia/Bindings/CMakeLists.txt index 63714b8f3..81f09673b 100644 --- a/Kassiopeia/Bindings/CMakeLists.txt +++ b/Kassiopeia/Bindings/CMakeLists.txt @@ -118,6 +118,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 @@ -343,6 +344,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 From 3c38e14fe57726a1282ad642135df43d827ca0bb Mon Sep 17 00:00:00 2001 From: juniorpena <72262395+juniorpena@users.noreply.github.com> Date: Mon, 22 Aug 2022 17:06:53 -0400 Subject: [PATCH 07/22] Update KSIntCalculatorMott.h --- Kassiopeia/Interactions/Include/KSIntCalculatorMott.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Kassiopeia/Interactions/Include/KSIntCalculatorMott.h b/Kassiopeia/Interactions/Include/KSIntCalculatorMott.h index b92c59631..e90374f6e 100644 --- a/Kassiopeia/Interactions/Include/KSIntCalculatorMott.h +++ b/Kassiopeia/Interactions/Include/KSIntCalculatorMott.h @@ -9,7 +9,7 @@ * KSintCalculatorMott.h * * Date: August 22, 2022 - * Author: Junior Peña (juniorpe) + * Author: Junior Peña (juniorpena) */ namespace Kassiopeia From 9448d42dc1033f55b7ff8588aadd352cf2787249 Mon Sep 17 00:00:00 2001 From: juniorpena <72262395+juniorpena@users.noreply.github.com> Date: Mon, 6 May 2024 15:10:58 -0700 Subject: [PATCH 08/22] Update KSIntCalculatorMottBuilder.h --- .../Include/KSIntCalculatorMottBuilder.h | 34 +++++++++++++++++-- 1 file changed, 32 insertions(+), 2 deletions(-) diff --git a/Kassiopeia/Bindings/Interactions/Include/KSIntCalculatorMottBuilder.h b/Kassiopeia/Bindings/Interactions/Include/KSIntCalculatorMottBuilder.h index 6f6880f2c..a970146ba 100644 --- a/Kassiopeia/Bindings/Interactions/Include/KSIntCalculatorMottBuilder.h +++ b/Kassiopeia/Bindings/Interactions/Include/KSIntCalculatorMottBuilder.h @@ -18,11 +18,41 @@ template<> inline bool KSIntCalculatorMottBuilder::AddAttribute(KContainer* aCon } if (aContainer->GetName() == "theta_min") { aContainer->CopyTo(fObject, &KSIntCalculatorMott::SetThetaMin); - return true; + + 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); - return true; + 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") != 0) && (fObject->GetNucleus().compare("Ne") != 0)) { + initmsg(eError) << "\"" << fObject->GetNucleus() + << R"(" is not available for Mott scattering! Available gases: "He", "Ne")" << eom; + + return false; + } + else { + return true; + } } return false; } From b52ada7ce41d160402198f7b57524f9baa44ecf6 Mon Sep 17 00:00:00 2001 From: juniorpena <72262395+juniorpena@users.noreply.github.com> Date: Mon, 6 May 2024 15:11:36 -0700 Subject: [PATCH 09/22] Update KSIntCalculatorMottBuilder.cxx --- .../Interactions/Source/KSIntCalculatorMottBuilder.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Kassiopeia/Bindings/Interactions/Source/KSIntCalculatorMottBuilder.cxx b/Kassiopeia/Bindings/Interactions/Source/KSIntCalculatorMottBuilder.cxx index 7590801d0..58891e168 100644 --- a/Kassiopeia/Bindings/Interactions/Source/KSIntCalculatorMottBuilder.cxx +++ b/Kassiopeia/Bindings/Interactions/Source/KSIntCalculatorMottBuilder.cxx @@ -12,7 +12,8 @@ template<> KSIntCalculatorMottBuilder::~KComplexElement() = default; STATICINT sKSIntCalculatorMottStructure = KSIntCalculatorMottBuilder::Attribute("name") + KSIntCalculatorMottBuilder::Attribute("theta_min") + - KSIntCalculatorMottBuilder::Attribute("theta_max"); + KSIntCalculatorMottBuilder::Attribute("theta_max") + + KSIntCalculatorMottBuilder::Attribute("nucleus"); STATICINT sToolboxKSIntCalculatorMott = KSRootBuilder::ComplexElement("ksint_calculator_mott"); From 26034d4a481342f8254b9e31e8e271097d07bf0b Mon Sep 17 00:00:00 2001 From: juniorpena <72262395+juniorpena@users.noreply.github.com> Date: Mon, 6 May 2024 15:12:20 -0700 Subject: [PATCH 10/22] Update KSIntCalculatorMott.h --- .../Include/KSIntCalculatorMott.h | 70 +++++++++++++++---- 1 file changed, 55 insertions(+), 15 deletions(-) diff --git a/Kassiopeia/Interactions/Include/KSIntCalculatorMott.h b/Kassiopeia/Interactions/Include/KSIntCalculatorMott.h index e90374f6e..3bdd9c389 100644 --- a/Kassiopeia/Interactions/Include/KSIntCalculatorMott.h +++ b/Kassiopeia/Interactions/Include/KSIntCalculatorMott.h @@ -3,13 +3,12 @@ #include "KField.h" #include "KSIntCalculator.h" -#include "TF1.h" /* * KSintCalculatorMott.h * * Date: August 22, 2022 - * Author: Junior Peña (juniorpena) + * Author: Junior Ivan Peña (juniorpe) */ namespace Kassiopeia @@ -18,7 +17,7 @@ namespace Kassiopeia /* * The xml configuration for using this calculator is as follows: * - * where [lower_bound] is the lower limit, in radians, on the range of allowed scattering + * 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. */ @@ -40,11 +39,14 @@ class KSIntCalculatorMott : public KSComponentTemplate 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] * */ - void InitializeMDCS(const double E0); + double Normalized_RDCS(double theta); /** - * \brief Deinitializes Mott Differential Cross Section and Mott Total Cross Section + * \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] * */ - void DeinitializeMDCS(); + double Normalized_RTCSInverse(double x); - protected: - TF1* fMDCS; - TF1* fMTCS; }; } // namespace Kassiopeia From 341c15c8459c050dcbf0f49761a9be71fd6d54b1 Mon Sep 17 00:00:00 2001 From: juniorpena <72262395+juniorpena@users.noreply.github.com> Date: Mon, 6 May 2024 15:13:10 -0700 Subject: [PATCH 11/22] Update KSIntCalculatorMott.cxx --- .../Source/KSIntCalculatorMott.cxx | 178 +++++++++++------- 1 file changed, 107 insertions(+), 71 deletions(-) diff --git a/Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx b/Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx index c9d2b2553..9dfb77b6f 100644 --- a/Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx +++ b/Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx @@ -1,4 +1,6 @@ #include "KSIntCalculatorMott.h" +#include +#include #include "KRandom.h" using katrin::KRandom; @@ -10,13 +12,12 @@ using KGeoBag::KThreeVector; namespace Kassiopeia { -KSIntCalculatorMott::KSIntCalculatorMott() : fThetaMin(0.), fThetaMax(0.), fMDCS(nullptr), fMTCS(nullptr) {} +KSIntCalculatorMott::KSIntCalculatorMott() : fThetaMin(0.), fThetaMax(0.), fNucleus("") {} KSIntCalculatorMott::KSIntCalculatorMott(const KSIntCalculatorMott& aCopy) : KSComponent(aCopy), fThetaMin(aCopy.fThetaMin), fThetaMax(aCopy.fThetaMax), - fMDCS(nullptr), - fMTCS(nullptr) + fNucleus(aCopy.fNucleus) {} KSIntCalculatorMott* KSIntCalculatorMott::Clone() const { @@ -28,14 +29,10 @@ void KSIntCalculatorMott::CalculateCrossSection(const KSParticle& anInitialParti { double tInitialKineticEnergy = anInitialParticle.GetKineticEnergy_eV(); - InitializeMDCS(tInitialKineticEnergy); - double fCrossSection = (fMTCS->Eval(fThetaMax) - fMTCS->Eval(fThetaMin)) * pow(10, -30); - - DeinitializeMDCS(); + double fCrossSection = (MTCS(fThetaMax, tInitialKineticEnergy) - MTCS(fThetaMin, tInitialKineticEnergy)); aCrossSection = fCrossSection; - return; } void KSIntCalculatorMott::ExecuteInteraction(const KSParticle& anInitialParticle, KSParticle& aFinalParticle, @@ -51,7 +48,7 @@ void KSIntCalculatorMott::ExecuteInteraction(const KSParticle& anInitialParticle double tPhi = KRandom::GetInstance().Uniform(0., 2. * katrin::KConst::Pi()); - tMomentum.SetPolarAngle(tTheta); + tMomentum.SetPolarAngle(tTheta * katrin::KConst::Pi() / 180); tMomentum.SetAzimuthalAngle(tPhi); aFinalParticle.SetTime(tTime); @@ -60,7 +57,7 @@ void KSIntCalculatorMott::ExecuteInteraction(const KSParticle& anInitialParticle aFinalParticle.SetKineticEnergy_eV(tInitialKineticEnergy - tLostKineticEnergy); aFinalParticle.SetLabel(GetName()); - fStepAngularChange = tTheta * 180. / katrin::KConst::Pi(); + fStepAngularChange = tTheta; return; } @@ -69,11 +66,36 @@ double KSIntCalculatorMott::GetTheta(const double& anEnergy){ double tTheta; - InitializeMDCS(anEnergy); + std::vector theta_arr; - tTheta = fMDCS->GetRandom(fThetaMin, fThetaMax); + // Populate theta_arr + for (int i = 0; i <= 1000; ++i) { + theta_arr.push_back(fThetaMin + i * (fThetaMax - fThetaMin) / 1000.0); + } - DeinitializeMDCS(); + double k = 0.0; + + // Calculating rescaling factor for rejection sampling + for (double theta : theta_arr) { + double MDCS_val = MDCS(theta, anEnergy); + double PDF_val = Normalized_RDCS(theta); + double ratio = MDCS_val / PDF_val; + if (ratio > k) { + k = MDCS_val / PDF_val; + } else { + break; + } + } + + // 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; @@ -82,94 +104,108 @@ double KSIntCalculatorMott::GetTheta(const double& anEnergy){ double KSIntCalculatorMott::GetEnergyLoss(const double& anEnergy, const double& aTheta){ double M, me, p, anELoss; - M = 4.0026 * katrin::KConst::AtomicMassUnit_eV(); // eV/c^2 + double amu = 0.0; + + if (fNucleus == "He") { + amu = 4.0026; + } else if (fNucleus == "Ne") { + amu = 20.1797; + } + + 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)); + 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; } -void KSIntCalculatorMott::InitializeMDCS(const double E0) -{ - - double beta = TMath::Sqrt( pow(E0, 2) + 2 * E0 * katrin::KConst::M_el_eV()) / (E0 + katrin::KConst::M_el_eV()); - /* Constants given in publication in header file. These are the constants for scattering off of Helium */ - static double b[5][6] = {{ 1. , 3.76476 * pow(10, -8), -3.05313 * pow(10, -7), -3.27422 * pow(10, -7), 2.44235 * pow(10, -6), 4.08754 * pow(10, -6)}, - { 2.35767 * pow(10, -2), 3.24642 * pow(10, -2), -6.37269 * pow(10, -4), -7.69160 * pow(10, -4), 5.28004 * pow(10, -3), 9.45642 * pow(10, -3)}, - { -2.73743 * pow(10, -1), -7.40767 * pow(10, -1), -4.98195 * pow(10, -1), 1.74337 * pow(10, -3), -1.25798 * pow(10, -2), -2.24046 * pow(10, -2)}, - { -7.79128 * pow(10, -4), -4.14495 * pow(10, -4), -1.62657 * pow(10, -3), -1.37286 * pow(10, -3), 1.04319 * pow(10, -2), 1.83488 * pow(10, -2)}, - { 2.02855 * pow(10, -4), 1.94598 * pow(10, -6), 4.30102 * pow(10, -4), 4.32180 * pow(10, -4), -3.31526 * pow(10, -3), -5.81788 * pow(10, -3)}}; - - double a[6] = {0, 0, 0, 0, 0, 0}; +double KSIntCalculatorMott::Beta(double const E0) { + return sqrt( pow(E0, 2) + 2 * E0 * katrin::KConst::M_el_eV()) / (E0 + katrin::KConst::M_el_eV()); +} - for (int j = 0; j < 5; j++){ - for (int k = 0; k < 6; k++){ - a[j] += b[j][k] * pow((beta - 0.7181287), k); +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 (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 (fNucleus == "Ne") { + a[5] = 10; + 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); } } - - /* ROOT TF1 takes in analytical formulas as strings. This is the conversion into strings */ - std::string Rmott; - std::ostringstream RmottStream; - - RmottStream << std::fixed; - RmottStream << std::setprecision(12); - RmottStream << a[0] << " + " << a[1] << " * pow( ( 1 - cos(x) ) , 0.5) + " << a[2] << " * ( 1 - cos(x) ) + " << a[3] << " * pow( ( 1 - cos(x) ) , 1.5) + " << a[4] << " * pow( ( 1 - cos(x) ) , 2)"; - Rmott = RmottStream.str(); + return a; +} +double KSIntCalculatorMott::MDCS(double theta, const double E0) { + double r_e = 2.8179403227 * pow(10, -15); - std::string RutherfordDCS; - std::ostringstream RutherfordDCSStream; + std::vector a = RMott_coeffs(E0); - RutherfordDCSStream << std::fixed; - RutherfordDCSStream << std::setprecision(12); + double Z = a[5]; - RutherfordDCSStream << "pow((2 * 2.8179403227 ), 2) * ( (1 - pow(" << beta << ", 2) ) / ( pow(" << beta << ", 4) )) * (1 / pow( (1 - cos(x)), 2))"; - RutherfordDCS = RutherfordDCSStream.str(); + return pow(Z, 2) * pow(r_e, 2) * ( (1 - pow(Beta(E0), 2))/(pow(Beta(E0), 4)) ) * ((a[0] * pow(1 - cos(theta), -2)) + (a[1] * pow((1 - cos(theta)), ((1 / 2) - 2))) + (a[2] * pow((1 - cos(theta)), (-1))) + (a[3] * pow((1 - cos(theta)), ((3 / 2) - 2))) + a[4] ); +} +double KSIntCalculatorMott::MTCS(double theta, const double E0) { + double r_e = 2.8179403227 * pow(10, -15); - std::string MottDCS = "(" + RutherfordDCS + ") * (" + Rmott + ")"; + std::vector a = RMott_coeffs(E0); + double Z = a[5]; - std::string MottTCS; - std::ostringstream MottTCSStream; - std::ostringstream MottTCSStreamThetaIndepPrefactor; - MottTCSStream << std::fixed; - MottTCSStream << std::setprecision(12); - - MottTCSStreamThetaIndepPrefactor << std::fixed; - MottTCSStreamThetaIndepPrefactor << std::setprecision(12); + 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))); +} - MottTCSStreamThetaIndepPrefactor << "pow((2 * 2.8179403227 ), 2) * ( (1 - pow(" << beta << ", 2) ) / ( pow(" << beta << ", 4) )) * "; +double KSIntCalculatorMott::Normalized_RDCS(double theta) { - MottTCSStream << "(2 * " << katrin::KConst::Pi() << " * (-4 * sqrt(2) * ( (" << a[1] << ") + (-1) * (" << a[2] - << ") * sqrt(( 1 - cos(x) )) * log( sin(x / 2) )) * pow(sin(x / 2), 4) + 4 * sqrt(2) * (2 * (" - << a[3] << ") + (" << a[4] << ") * sqrt( (1 - cos(x))) ) * pow(sin(x / 2), 6) - 2 * (" << a[0] - << ") * pow(sin( x / 2 ), 3) ) ) / ( pow((-1 + cos(x)), 2) * sin( x / 2 ))"; + double C = 1 / ((1 / (1 - cos(fThetaMin* katrin::KConst::Pi() / 180))) - (1 / (1 - cos(fThetaMax* katrin::KConst::Pi() / 180)))); - MottTCS = MottTCSStreamThetaIndepPrefactor.str() + MottTCSStream.str(); + return C / pow((1 - cos(theta* katrin::KConst::Pi() / 180)), 2); +} - fMDCS = new TF1("function", MottDCS.c_str(), fThetaMin, fThetaMax); - fMTCS = new TF1("function", MottTCS.c_str(), fThetaMin, fThetaMax); +double KSIntCalculatorMott::Normalized_RTCSInverse(double x) { - return; -} + double C = 1 / ((1 / (1 - cos(fThetaMin * katrin::KConst::Pi() / 180))) - (1 / (1 - cos(fThetaMax * katrin::KConst::Pi() / 180)))); -void KSIntCalculatorMott::DeinitializeMDCS() -{ - delete fMDCS; - fMDCS = nullptr; - delete fMTCS; - fMTCS = nullptr; - return; + return acos( 1 - (1 / ((1 / (1 - cos(fThetaMin * katrin::KConst::Pi() / 180))) - x / C)) ) * 180 / katrin::KConst::Pi(); } } // namespace Kassiopeia From 335814e0698c115a5294c9a081351126bd4ae4ff Mon Sep 17 00:00:00 2001 From: juniorpena <72262395+juniorpena@users.noreply.github.com> Date: Tue, 7 May 2024 11:18:54 -0700 Subject: [PATCH 12/22] Update Interactions.xml --- Kassiopeia/XML/Complete/Interactions.xml | 29 ++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/Kassiopeia/XML/Complete/Interactions.xml b/Kassiopeia/XML/Complete/Interactions.xml index 0e92da518..d5b66d6d1 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. --> + + + + + + + From 08659037cb1e704fadbc32853043210bba64213f Mon Sep 17 00:00:00 2001 From: juniorpena <72262395+juniorpena@users.noreply.github.com> Date: Mon, 10 Jun 2024 11:23:26 -0700 Subject: [PATCH 13/22] Update Kassiopeia/Bindings/Interactions/Include/KSIntCalculatorMottBuilder.h Co-authored-by: 2xB <31772910+2xB@users.noreply.github.com> --- .../Bindings/Interactions/Include/KSIntCalculatorMottBuilder.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Kassiopeia/Bindings/Interactions/Include/KSIntCalculatorMottBuilder.h b/Kassiopeia/Bindings/Interactions/Include/KSIntCalculatorMottBuilder.h index a970146ba..d9738fb9c 100644 --- a/Kassiopeia/Bindings/Interactions/Include/KSIntCalculatorMottBuilder.h +++ b/Kassiopeia/Bindings/Interactions/Include/KSIntCalculatorMottBuilder.h @@ -19,7 +19,7 @@ template<> inline bool KSIntCalculatorMottBuilder::AddAttribute(KContainer* aCon if (aContainer->GetName() == "theta_min") { aContainer->CopyTo(fObject, &KSIntCalculatorMott::SetThetaMin); - if ((fObject->GetThetaMin() == 0.0)) { + if ((fObject->GetThetaMin() <= 0.0)) { initmsg(eError) << "\"" << fObject->GetThetaMin() << R"(" is an invalid lower bound for Mott scattering! Change to a value > 0)" << eom; From 317fcae9c2cea9ef561a8684c493bfcb8b0d02b7 Mon Sep 17 00:00:00 2001 From: juniorpena <72262395+juniorpena@users.noreply.github.com> Date: Mon, 10 Jun 2024 11:33:55 -0700 Subject: [PATCH 14/22] Update KSIntCalculatorMott.cxx I've added code and comments for the choice of size of theta_arr responding to previous comment. I've also added comments wherever there are hard coded values, for where those values come from --- .../Source/KSIntCalculatorMott.cxx | 56 +++++++++++-------- 1 file changed, 34 insertions(+), 22 deletions(-) diff --git a/Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx b/Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx index 9dfb77b6f..bbe96e160 100644 --- a/Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx +++ b/Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx @@ -1,6 +1,6 @@ #include "KSIntCalculatorMott.h" #include -#include +#include #include "KRandom.h" using katrin::KRandom; @@ -51,12 +51,22 @@ void KSIntCalculatorMott::ExecuteInteraction(const KSParticle& anInitialParticle tMomentum.SetPolarAngle(tTheta * katrin::KConst::Pi() / 180); tMomentum.SetAzimuthalAngle(tPhi); + 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(tMomentum); + aFinalParticle.SetMomentum(tFinalDirection); aFinalParticle.SetKineticEnergy_eV(tInitialKineticEnergy - tLostKineticEnergy); aFinalParticle.SetLabel(GetName()); + fStepNInteractions = 1; + fStepEnergyLoss = tLostKineticEnergy; fStepAngularChange = tTheta; return; @@ -65,28 +75,30 @@ void KSIntCalculatorMott::ExecuteInteraction(const KSParticle& anInitialParticle 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 - std::vector theta_arr; - // Populate theta_arr - for (int i = 0; i <= 1000; ++i) { - theta_arr.push_back(fThetaMin + i * (fThetaMax - fThetaMin) / 1000.0); - } + resolution = ceil((fThetaMax - fThetaMin) * resolution_factor); + - double k = 0.0; + // 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) { - double MDCS_val = MDCS(theta, anEnergy); - double PDF_val = Normalized_RDCS(theta); - double ratio = MDCS_val / PDF_val; - if (ratio > k) { - k = MDCS_val / PDF_val; - } else { - break; - } + 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)); @@ -107,9 +119,9 @@ double KSIntCalculatorMott::GetEnergyLoss(const double& anEnergy, const double& double amu = 0.0; if (fNucleus == "He") { - amu = 4.0026; + amu = 4.0026; // mass in amu from W.J. Huang et al 2021 Chinese Phys. C 45 030002 } else if (fNucleus == "Ne") { - amu = 20.1797; + 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 @@ -140,7 +152,7 @@ std::vector KSIntCalculatorMott::RMott_coeffs(double const E0) { {-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 (fNucleus == "Ne") { - a[5] = 10; + 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 }, @@ -160,18 +172,18 @@ std::vector KSIntCalculatorMott::RMott_coeffs(double const E0) { double KSIntCalculatorMott::MDCS(double theta, const double E0) { - double r_e = 2.8179403227 * pow(10, -15); + 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), -2)) + (a[1] * pow((1 - cos(theta)), ((1 / 2) - 2))) + (a[2] * pow((1 - cos(theta)), (-1))) + (a[3] * pow((1 - cos(theta)), ((3 / 2) - 2))) + a[4] ); + 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); + double r_e = 2.8179403227 * pow(10, -15); // classical electron radius std::vector a = RMott_coeffs(E0); From 811216646e7e371cccdb4a89e89ed4d94d01c7c2 Mon Sep 17 00:00:00 2001 From: juniorpena <72262395+juniorpena@users.noreply.github.com> Date: Tue, 16 Jul 2024 11:36:12 -0700 Subject: [PATCH 15/22] Update KSIntCalculatorMott.cxx Merged main branch, and fixed bug causing Ubuntu 20.04 CI check failure --- Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx b/Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx index bbe96e160..1ee8dfc92 100644 --- a/Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx +++ b/Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx @@ -7,7 +7,8 @@ using katrin::KRandom; #include "KConst.h" -using KGeoBag::KThreeVector; +#include "KThreeVector.hh" +using katrin::KThreeVector; namespace Kassiopeia { From 9d9f6689e39c25c21d1eaa1c71cc9315b4796585 Mon Sep 17 00:00:00 2001 From: juniorpena <72262395+juniorpena@users.noreply.github.com> Date: Thu, 18 Jul 2024 12:39:29 -0700 Subject: [PATCH 16/22] Update KSIntCalculatorMott.cxx Found a bug after running simulations, the momentum of the final state particle was being changed twice when executing the interaction --- Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx | 2 -- 1 file changed, 2 deletions(-) diff --git a/Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx b/Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx index 1ee8dfc92..72beece50 100644 --- a/Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx +++ b/Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx @@ -49,8 +49,6 @@ void KSIntCalculatorMott::ExecuteInteraction(const KSParticle& anInitialParticle double tPhi = KRandom::GetInstance().Uniform(0., 2. * katrin::KConst::Pi()); - tMomentum.SetPolarAngle(tTheta * katrin::KConst::Pi() / 180); - tMomentum.SetAzimuthalAngle(tPhi); KThreeVector tOrthogonalOne = tMomentum.Orthogonal(); KThreeVector tOrthogonalTwo = tMomentum.Cross(tOrthogonalOne); From c66c2f66a7cbfcd73247902f5e0184787c1737ae Mon Sep 17 00:00:00 2001 From: juniorpena <72262395+juniorpena@users.noreply.github.com> Date: Thu, 12 Sep 2024 13:11:33 -0700 Subject: [PATCH 17/22] Update KSIntCalculatorMott.cxx --- Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx b/Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx index 72beece50..a93eeae71 100644 --- a/Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx +++ b/Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx @@ -117,9 +117,9 @@ double KSIntCalculatorMott::GetEnergyLoss(const double& anEnergy, const double& double amu = 0.0; - if (fNucleus == "He") { + 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") { + } else if (fNucleus == "Ne-20") { amu = 19.9924; // mass in amu from W.J. Huang et al 2021 Chinese Phys. C 45 030002 } @@ -143,14 +143,14 @@ 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 (fNucleus == "He") { + if (fNucleus == "He-4") { // independent of mass number, dependent on atomic number (ie. He-3 would have same data as He-4) 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 (fNucleus == "Ne") { + } else if (fNucleus == "Ne-20") { // independent of mass number, dependent on atomic number (ie. Ne-21 would have same data as Ne-20) 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}, From 58d7ac0a5e25dd5887367c3a9f91d39ded236452 Mon Sep 17 00:00:00 2001 From: juniorpena <72262395+juniorpena@users.noreply.github.com> Date: Thu, 12 Sep 2024 13:12:15 -0700 Subject: [PATCH 18/22] Update KSIntCalculatorMottBuilder.h --- .../Interactions/Include/KSIntCalculatorMottBuilder.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Kassiopeia/Bindings/Interactions/Include/KSIntCalculatorMottBuilder.h b/Kassiopeia/Bindings/Interactions/Include/KSIntCalculatorMottBuilder.h index d9738fb9c..e843f285c 100644 --- a/Kassiopeia/Bindings/Interactions/Include/KSIntCalculatorMottBuilder.h +++ b/Kassiopeia/Bindings/Interactions/Include/KSIntCalculatorMottBuilder.h @@ -44,9 +44,9 @@ template<> inline bool KSIntCalculatorMottBuilder::AddAttribute(KContainer* aCon if (aContainer->GetName() == "nucleus") { aContainer->CopyTo(fObject, &KSIntCalculatorMott::SetNucleus); - if ((fObject->GetNucleus().compare("He") != 0) && (fObject->GetNucleus().compare("Ne") != 0)) { + 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", "Ne")" << eom; + << R"(" is not available for Mott scattering! Available gases: "He-4", "Ne-20")" << eom; return false; } From c338bf90eaac908defcca1b538db3e68f15db062 Mon Sep 17 00:00:00 2001 From: juniorpena <72262395+juniorpena@users.noreply.github.com> Date: Thu, 12 Sep 2024 13:15:44 -0700 Subject: [PATCH 19/22] Update Kassiopeia/XML/Complete/Interactions.xml Co-authored-by: 2xB <31772910+2xB@users.noreply.github.com> --- Kassiopeia/XML/Complete/Interactions.xml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Kassiopeia/XML/Complete/Interactions.xml b/Kassiopeia/XML/Complete/Interactions.xml index d5b66d6d1..5cd7cd4aa 100644 --- a/Kassiopeia/XML/Complete/Interactions.xml +++ b/Kassiopeia/XML/Complete/Interactions.xml @@ -114,11 +114,11 @@ - +