Code indexing in gitaly is broken and leads to code not being visible to the user. We work on the issue with highest priority.

Skip to content
Snippets Groups Projects

Resolve "Reviewing physics behind particle matter interaction models"

Merged Resolve "Reviewing physics behind particle matter interaction models"
634-reviewing-physics-behind-particle-matter-interaction-models into master
Merged ext-calvo_p requested to merge 634-reviewing-physics-behind-particle-matter-interaction-models into master
Compare
and
Show latest version
2 files
+ 63
− 65
Compare changes
  • Side-by-side
  • Inline
Files
2
src/Classic/Solvers/CollimatorPhysics.cpp
+ 48
− 50
@@ -48,49 +48,49 @@
namespace {
struct DegraderInsideTester: public InsideTester {
explicit DegraderInsideTester(ElementBase * el) {
explicit DegraderInsideTester(ElementBase* el) {
deg_m = static_cast<Degrader*>(el);
}
virtual
bool checkHit(const Vector_t &R) override {
bool checkHit(const Vector_t& R) override {
return deg_m->isInMaterial(R(2));
}
private:
Degrader *deg_m;
Degrader* deg_m;
};
struct CollimatorInsideTester: public InsideTester {
explicit CollimatorInsideTester(ElementBase * el) {
explicit CollimatorInsideTester(ElementBase* el) {
col_m = static_cast<CCollimator*>(el);
}
virtual
bool checkHit(const Vector_t &R) override {
bool checkHit(const Vector_t& R) override {
return col_m->checkPoint(R(0), R(1));
}
private:
CCollimator *col_m;
CCollimator* col_m;
};
struct FlexCollimatorInsideTester: public InsideTester {
explicit FlexCollimatorInsideTester(ElementBase * el) {
explicit FlexCollimatorInsideTester(ElementBase* el) {
col_m = static_cast<FlexibleCollimator*>(el);
}
virtual
bool checkHit(const Vector_t &R) override {
bool checkHit(const Vector_t& R) override {
return col_m->isStopped(R);
}
private:
FlexibleCollimator *col_m;
FlexibleCollimator* col_m;
};
}
CollimatorPhysics::CollimatorPhysics(const std::string &name,
ElementBase *element,
std::string &material,
CollimatorPhysics::CollimatorPhysics(const std::string& name,
ElementBase* element,
std::string& material,
bool enableRutherford,
double lowEnergyThr):
ParticleMatterInteractionHandler(name, element),
@@ -183,8 +183,8 @@ void CollimatorPhysics::configureMaterialParameters() {
A5_c = material->getStoppingPowerFitCoefficients(Physics::Material::A5);
}
void CollimatorPhysics::apply(PartBunchBase<double, 3> *bunch,
const std::pair<Vector_t, double> &boundingSphere) {
void CollimatorPhysics::apply(PartBunchBase<double, 3>* bunch,
const std::pair<Vector_t, double>& boundingSphere) {
IpplTimings::startTimer(DegraderApplyTimer_m);
/*
@@ -257,9 +257,9 @@ void CollimatorPhysics::computeInteraction() {
for (size_t i = 0; i < locParts_m.size(); ++i) {
if (locParts_m[i].label != -1) {
Vector_t &R = locParts_m[i].Rincol;
Vector_t &P = locParts_m[i].Pincol;
double &dt = locParts_m[i].DTincol;
Vector_t& R = locParts_m[i].Rincol;
Vector_t& P = locParts_m[i].Pincol;
double& dt = locParts_m[i].DTincol;
if (hitTester_m->checkHit(R)) {
bool pdead = computeEnergyLoss(P, dt);
@@ -293,7 +293,7 @@ void CollimatorPhysics::computeInteraction() {
/// See Particle Physics Booklet, chapter 'Passage of particles through matter' or
/// Review of Particle Physics, DOI: 10.1103/PhysRevD.86.010001, page 329 ff
// -------------------------------------------------------------------------
bool CollimatorPhysics::computeEnergyLoss(Vector_t &P,
bool CollimatorPhysics::computeEnergyLoss(Vector_t& P,
const double deltat,
bool includeFluctuations) const {
@@ -318,7 +318,7 @@ bool CollimatorPhysics::computeEnergyLoss(Vector_t &P,
const double deltasrho = deltas * m2cm * rho_m;
if (Ekin >= 600) {
constexpr double massRatio = massElectron_keV / mass_keV;
const double massRatio = massElectron_keV / mass_keV;
double Tmax = (2.0 * massElectron_keV * betaSqr * gammaSqr /
(std::pow(gamma + massRatio, 2) - (gammaSqr - 1.0)));
@@ -329,8 +329,8 @@ bool CollimatorPhysics::computeEnergyLoss(Vector_t &P,
} else if (Ekin > 10) {
constexpr double massProton_amu = Physics::m_p / Physics::amu;
const double Ts = Ekin / massProton_amu;
const double epsilon_low = A2_c * pow(Ts, 0.45);
const double epsilon_high = (A3_c / Ts) * log(1 + (A4_c / Ts) + (A5_c * Ts));
const double epsilon_low = A2_c * std::pow(Ts, 0.45);
const double epsilon_high = (A3_c / Ts) * std::log(1 + (A4_c / Ts) + (A5_c * Ts));
const double epsilon = (epsilon_low * epsilon_high) / (epsilon_low + epsilon_high);
dEdx = -epsilon / (1e18 * (A_m / Physics::Avo));
@@ -355,35 +355,35 @@ bool CollimatorPhysics::computeEnergyLoss(Vector_t &P,
// For details see:
// J. Beringer et al. (Particle Data Group), "Passage of particles through matter"
// Phys. Rev. D 86, 010001 (2012),
void CollimatorPhysics::applyRotation(Vector_t &P,
Vector_t &R,
void CollimatorPhysics::applyRotation(Vector_t& P,
Vector_t& R,
double shift,
double thetacou) {
// Calculate the angle between the transverse and longitudinal component of the momentum
double Psixz = std::fmod(std::atan2(P(0), P(2)) + Physics::two_pi, Physics::two_pi);
R(0) = R(0) + shift * cos(Psixz);
R(2) = R(2) - shift * sin(Psixz);
R(0) = R(0) + shift * std::cos(Psixz);
R(2) = R(2) - shift * std::sin(Psixz);
// Apply the rotation about the random angle thetacou
double Px = P(0);
P(0) = Px * cos(thetacou) + P(2) * sin(thetacou);
P(2) = -Px * sin(thetacou) + P(2) * cos(thetacou);
P(0) = Px * std::cos(thetacou) + P(2) * std::sin(thetacou);
P(2) = -Px * std::sin(thetacou) + P(2) * std::cos(thetacou);
}
void CollimatorPhysics::applyRandomRotation(Vector_t &P, double theta0) {
void CollimatorPhysics::applyRandomRotation(Vector_t& P, double theta0) {
double thetaru = 2.5 / std::sqrt(gsl_rng_uniform(rGen_m)) * 2.0 * theta0;
double phiru = Physics::two_pi * gsl_rng_uniform(rGen_m);
double normPtrans = std::sqrt(P(0) * P(0) + P(1) * P(1));
double Theta = std::atan(normPtrans / std::abs(P(2)));
double CosT = cos(Theta);
double SinT = sin(Theta);
double CosT = std::cos(Theta);
double SinT = std::sin(Theta);
Vector_t X(cos(phiru)*sin(thetaru),
sin(phiru)*sin(thetaru),
cos(thetaru));
Vector_t X(std::cos(phiru)*std::sin(thetaru),
std::sin(phiru)*std::sin(thetaru),
std::cos(thetaru));
X *= euclidean_norm(P);
Vector_t W(-P(1), P(0), 0.0);
@@ -396,11 +396,10 @@ void CollimatorPhysics::applyRandomRotation(Vector_t &P, double theta0) {
/// Coulomb Scattering: Including Multiple Coulomb Scattering and large angle Rutherford Scattering.
/// Using the distribution given in Classical Electrodynamics, by J. D. Jackson.
//--------------------------------------------------------------------------
void CollimatorPhysics::computeCoulombScattering(Vector_t &R,
Vector_t &P,
void CollimatorPhysics::computeCoulombScattering(Vector_t& R,
Vector_t& P,
double dt) {
constexpr double GeV2eV = 1e9;
constexpr double sqrtThreeInv = 0.57735026918962576451; // sqrt(1.0 / 3.0)
const double normP = euclidean_norm(P);
const double gammaSqr = std::pow(normP, 2) + 1.0;
@@ -408,7 +407,7 @@ void CollimatorPhysics::computeCoulombScattering(Vector_t &R,
const double deltas = dt * beta * Physics::c;
const double theta0 = (13.6e6 / (beta * normP * mass_m) *
charge_m * std::sqrt(deltas / X0_m) *
(1.0 + 0.038 * log(deltas / X0_m)));
(1.0 + 0.038 * std::log(deltas / X0_m)));
double phi = Physics::two_pi * gsl_rng_uniform(rGen_m);
for (unsigned int i = 0; i < 2; ++ i) {
@@ -441,7 +440,7 @@ void CollimatorPhysics::computeCoulombScattering(Vector_t &R,
}
}
void CollimatorPhysics::addBackToBunch(PartBunchBase<double, 3> *bunch) {
void CollimatorPhysics::addBackToBunch(PartBunchBase<double, 3>* bunch) {
const size_t nL = locParts_m.size();
if (nL == 0) return;
@@ -450,7 +449,7 @@ void CollimatorPhysics::addBackToBunch(PartBunchBase<double, 3> *bunch) {
const double elementLength = element_ref_m->getElementLength();
for (size_t i = 0; i < nL; ++ i) {
Vector_t &R = locParts_m[i].Rincol;
Vector_t& R = locParts_m[i].Rincol;
if (R[2] >= elementLength) {
@@ -460,7 +459,6 @@ void CollimatorPhysics::addBackToBunch(PartBunchBase<double, 3> *bunch) {
Binincol is still <0, but now the particle is rediffused
from the material and hence this is not a "lost" particle anymore
*/
bunch->Bin[numLocalParticles] = 1;
bunch->R[numLocalParticles] = R;
bunch->P[numLocalParticles] = locParts_m[i].Pincol;
@@ -487,8 +485,8 @@ void CollimatorPhysics::addBackToBunch(PartBunchBase<double, 3> *bunch) {
deleteParticleFromLocalVector();
}
void CollimatorPhysics::copyFromBunch(PartBunchBase<double, 3> *bunch,
const std::pair<Vector_t, double> &boundingSphere)
void CollimatorPhysics::copyFromBunch(PartBunchBase<double, 3>* bunch,
const std::pair<Vector_t, double>& boundingSphere)
{
const size_t nL = bunch->getLocalNum();
if (nL == 0) return;
@@ -558,8 +556,8 @@ void CollimatorPhysics::print(Inform &msg) {
Inform::FmtFlags_t ff = msg.flags();
if (totalPartsInMat_m > 0 ||
bunchToMatStat_m > 0 ||
rediffusedStat_m > 0 ||
bunchToMatStat_m > 0 ||
rediffusedStat_m > 0 ||
stoppedPartStat_m > 0) {
OPALTimer::Timer time;
@@ -614,8 +612,8 @@ void CollimatorPhysics::deleteParticleFromLocalVector() {
void CollimatorPhysics::push() {
for (size_t i = 0; i < locParts_m.size(); ++i) {
Vector_t &R = locParts_m[i].Rincol;
Vector_t &P = locParts_m[i].Pincol;
Vector_t& R = locParts_m[i].Rincol;
Vector_t& P = locParts_m[i].Pincol;
double gamma = Util::getGamma(P);
R += 0.5 * dT_m * Physics::c * P / gamma;
@@ -632,9 +630,9 @@ void CollimatorPhysics::setTimeStepForLeavingParticles() {
const double elementLength = element_ref_m->getElementLength();
for (size_t i = 0; i < locParts_m.size(); ++i) {
Vector_t &R = locParts_m[i].Rincol;
Vector_t &P = locParts_m[i].Pincol;
double &dt = locParts_m[i].DTincol;
Vector_t& R = locParts_m[i].Rincol;
Vector_t& P = locParts_m[i].Pincol;
double& dt = locParts_m[i].DTincol;
double gamma = Util::getGamma(P);
Vector_t stepLength = dT_m * Physics::c * P / gamma;
@@ -654,7 +652,7 @@ void CollimatorPhysics::setTimeStepForLeavingParticles() {
void CollimatorPhysics::resetTimeStep() {
for (size_t i = 0; i < locParts_m.size(); ++i) {
double &dt = locParts_m[i].DTincol;
double& dt = locParts_m[i].DTincol;
dt = dT_m;
}