From 199853018a76aa8c28c786279a2d9eace0e0339c Mon Sep 17 00:00:00 2001 From: muffinfan Date: Wed, 24 Jan 2024 10:41:59 +0100 Subject: [PATCH 1/5] fix sampling from singly differential cross section in hydrogen ionization --- .../Source/KSIntCalculatorHydrogen.cxx | 68 +++++++++++++------ 1 file changed, 48 insertions(+), 20 deletions(-) diff --git a/Kassiopeia/Interactions/Source/KSIntCalculatorHydrogen.cxx b/Kassiopeia/Interactions/Source/KSIntCalculatorHydrogen.cxx index 47b843262..0d4791502 100644 --- a/Kassiopeia/Interactions/Source/KSIntCalculatorHydrogen.cxx +++ b/Kassiopeia/Interactions/Source/KSIntCalculatorHydrogen.cxx @@ -1451,33 +1451,61 @@ void KSIntCalculatorHydrogenIonisation::ExecuteInteraction(const KSParticle& anI void KSIntCalculatorHydrogenIonisation::CalculateFinalEnergy(const double aReducedInitalEnergy, double& aReducedFinalEnergy) { - double tReducedFinalEnergy; - double tCrossSection = 0; - - - double sigma_max; - CalculateEnergyDifferentialCrossSection(aReducedInitalEnergy, aReducedInitalEnergy - 1.01, sigma_max); - double sigma_min; - CalculateEnergyDifferentialCrossSection(aReducedInitalEnergy, (aReducedInitalEnergy - 1) / 2, sigma_min); + // Get final energy of primary electron after ionization scattering via + // rejection sampling from the SDCS with a custom powerlaw proposal distribution: + // pdf(e) = 2*sigma_min*1/(e+1)^2.5, where sigma_min = SCDS(eMin) + // (compatible with shape of Rudd 1991 SDCS parametrization). + // Proposal samples are generated through inverse transform sampling. + + // local output variable + double tReducedFinalEnergy; + // relevant energy values + //double e_max = (aReducedInitalEnergy-(1+1e-12)); + double e_hlf = (aReducedInitalEnergy-1.)/2.; + //double e_min = (0+1e-12); + + // value of SDCS at those values + //double sigma_max; CalculateEnergyDifferentialCrossSection(aReducedInitalEnergy, aReducedInitalEnergy-(1.+1e-12), sigma_max); + //double sigma_hlf; CalculateEnergyDifferentialCrossSection(aReducedInitalEnergy, (aReducedInitalEnergy-1.)/2., sigma_hlf); + double sigma_min; CalculateEnergyDifferentialCrossSection(aReducedInitalEnergy, 0.+1e-12, sigma_min); + + // normalization factor for ppf + double norm = (2./3.*(1-1/pow((e_hlf+1), 3./2.))); + + // pdf of proposal distribution defined on e in [e_min,e_hlf] + auto _pdf = [=] (double e) { return 2*sigma_min*1/pow((e+1.),2.5); }; + + // ppf of proposal distribution defined on u in [0,1] + auto _ppf = [=] (double u) { return (pow(2.,(2./3.))-pow((2.-3.*u*norm),(2./3.)))/(pow((2.-3.*u*norm),(2./3.))); }; + + // generate only energies below minimum at e_hlf, i.e. energy of secondary + // later exploit symmetry around e_hlf to get energy of primary electron + // maximal u value (ppf(u_max) = e_hlf) + double u_max = (2./3.)*(1-1/pow((e_hlf+1),(3./2.)))/norm; + + // variable that stores cross section at tried sample + double tCrossSection = 0.; + while (true) { - tReducedFinalEnergy = - KRandom::GetInstance().Uniform((aReducedInitalEnergy - 1.) / 2., aReducedInitalEnergy - 1., false, true); - + // sample energy from proposal distribution + double tUniform = KRandom::GetInstance().Uniform(0., u_max, false, true); + tReducedFinalEnergy = _ppf(tUniform); + + // get uniform random number between 0 and proposal pdf at sample energy + double tRandom = KRandom::GetInstance().Uniform(0., _pdf(tReducedFinalEnergy), false, true); + + // get value of target pdf for sample energy CalculateEnergyDifferentialCrossSection(aReducedInitalEnergy, tReducedFinalEnergy, tCrossSection); - - double tRandom = KRandom::GetInstance().Uniform(0., - 2. * (sigma_max - sigma_min) / (aReducedInitalEnergy - 1.) * - (tReducedFinalEnergy - (aReducedInitalEnergy - 1.) / 2.), - false, - true); - + + // accept or reject if (tRandom < tCrossSection) break; } - - aReducedFinalEnergy = tReducedFinalEnergy; + + // convert to energy of primary by subtracting binding energy and secondary energy from input + aReducedFinalEnergy = aReducedInitalEnergy-1.-tReducedFinalEnergy; } void KSIntCalculatorHydrogenIonisation::CalculateTheta(const double aReducedInitialEnergy, From daa5927d7d009cf281d771751aaf576a069c6b7e Mon Sep 17 00:00:00 2001 From: Richard Salomon Date: Mon, 11 Mar 2024 11:15:00 +0100 Subject: [PATCH 2/5] Reset KSTermOutput terminator value when terminator is called --- Kassiopeia/Terminators/Include/KSTermOutput.h | 1 + 1 file changed, 1 insertion(+) diff --git a/Kassiopeia/Terminators/Include/KSTermOutput.h b/Kassiopeia/Terminators/Include/KSTermOutput.h index 0e7d76383..e3961e918 100644 --- a/Kassiopeia/Terminators/Include/KSTermOutput.h +++ b/Kassiopeia/Terminators/Include/KSTermOutput.h @@ -42,6 +42,7 @@ template class KSTermOutput : public KSComponentTemplate= fMaxValue || *fValue <= fMinValue) { + *fValue = 0; aFlag = true; return; } From bb52aec9e4873605b673381451f3b557906f6e0b Mon Sep 17 00:00:00 2001 From: Martin Descher Date: Mon, 7 Oct 2024 14:02:41 +0200 Subject: [PATCH 3/5] ionization secondary angle fix by Julanan --- .../Interactions/Source/KSIntCalculatorHydrogen.cxx | 12 +++--------- 1 file changed, 3 insertions(+), 9 deletions(-) diff --git a/Kassiopeia/Interactions/Source/KSIntCalculatorHydrogen.cxx b/Kassiopeia/Interactions/Source/KSIntCalculatorHydrogen.cxx index 0d4791502..6c4b0c3b9 100644 --- a/Kassiopeia/Interactions/Source/KSIntCalculatorHydrogen.cxx +++ b/Kassiopeia/Interactions/Source/KSIntCalculatorHydrogen.cxx @@ -1427,15 +1427,9 @@ void KSIntCalculatorHydrogenIonisation::ExecuteInteraction(const KSParticle& anI fStepEnergyLoss = (tInitialKineticEnergy - tReducedFinalEnergy * BindingEnergy); // outgoing secondary - - tTheta = acos(KRandom::GetInstance().Uniform(-1., 1.)); - tPhi = KRandom::GetInstance().Uniform(0., 2. * katrin::KConst::Pi()); - - tOrthogonalOne = tInitialDirection.Orthogonal(); - tOrthogonalTwo = tInitialDirection.Cross(tOrthogonalOne); - tFinalDirection = tInitialDirection.Magnitude() * - (sin(tTheta) * (cos(tPhi) * tOrthogonalOne.Unit() + sin(tPhi) * tOrthogonalTwo.Unit()) + - cos(tTheta) * tInitialDirection.Unit()); + // energy: initial energy - primary energy - binding energy + // direction: use momentum conservation for free electron-electron scattering + tFinalDirection = tInitialDirection - tPrimaryDirection; KSParticle* tSecondary = KSParticleFactory::GetInstance().Create(11); (*tSecondary) = anInitialParticle; From 2340e6c42dc361b9f2b22fdd5601aa2f3e916ed7 Mon Sep 17 00:00:00 2001 From: Martin Descher Date: Tue, 8 Oct 2024 14:55:43 +0200 Subject: [PATCH 4/5] fix obvious mistake in previous commit, to be tested --- Kassiopeia/Interactions/Source/KSIntCalculatorHydrogen.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Kassiopeia/Interactions/Source/KSIntCalculatorHydrogen.cxx b/Kassiopeia/Interactions/Source/KSIntCalculatorHydrogen.cxx index 6c4b0c3b9..79b128360 100644 --- a/Kassiopeia/Interactions/Source/KSIntCalculatorHydrogen.cxx +++ b/Kassiopeia/Interactions/Source/KSIntCalculatorHydrogen.cxx @@ -1429,11 +1429,11 @@ void KSIntCalculatorHydrogenIonisation::ExecuteInteraction(const KSParticle& anI // outgoing secondary // energy: initial energy - primary energy - binding energy // direction: use momentum conservation for free electron-electron scattering - tFinalDirection = tInitialDirection - tPrimaryDirection; + KThreeVector tSecondaryDirection = tInitialDirection - aFinalParticle.GetMomentum(); KSParticle* tSecondary = KSParticleFactory::GetInstance().Create(11); (*tSecondary) = anInitialParticle; - tSecondary->SetMomentum(tFinalDirection); + tSecondary->SetMomentum(tSecondaryDirection); tSecondary->SetKineticEnergy_eV((tReducedInitialEnergy - tReducedFinalEnergy - 1.) * BindingEnergy); tSecondary->SetLabel(GetName()); From afe24150adb37726966db63adac070f7db19ed30 Mon Sep 17 00:00:00 2001 From: muffinfan Date: Thu, 12 Dec 2024 12:22:45 +0100 Subject: [PATCH 5/5] Update Kassiopeia/Terminators/Include/KSTermOutput.h Remove unrelated/unnecessary previous commit according to 2xB's suggestion. Co-authored-by: 2xB <31772910+2xB@users.noreply.github.com> --- Kassiopeia/Terminators/Include/KSTermOutput.h | 1 - 1 file changed, 1 deletion(-) diff --git a/Kassiopeia/Terminators/Include/KSTermOutput.h b/Kassiopeia/Terminators/Include/KSTermOutput.h index e3961e918..0e7d76383 100644 --- a/Kassiopeia/Terminators/Include/KSTermOutput.h +++ b/Kassiopeia/Terminators/Include/KSTermOutput.h @@ -42,7 +42,6 @@ template class KSTermOutput : public KSComponentTemplate= fMaxValue || *fValue <= fMinValue) { - *fValue = 0; aFlag = true; return; }