ScatteringPhysics.cpp 24.9 KB
Newer Older
1
//
2
// Class ScatteringPhysics
kraus's avatar
kraus committed
3
//   Defines the physical processes of beam scattering
4
//   and energy loss by heavy charged particles
gsell's avatar
gsell committed
5
//
6
// Copyright (c) 2009 - 2021, Bi, Yang, Stachel, Adelmann
7 8 9 10 11 12 13 14 15 16 17 18 19
//                            Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved.
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL.  If not, see <https://www.gnu.org/licenses/>.
//
kraus's avatar
kraus committed
20
#include "Solvers/ScatteringPhysics.h"
gsell's avatar
gsell committed
21
#include "Physics/Physics.h"
22
#include "Physics/Material.h"
23
#include "Algorithms/PartBunchBase.h"
24 25
#include "AbsBeamline/CCollimator.h"
#include "AbsBeamline/FlexibleCollimator.h"
adelmann's avatar
adelmann committed
26
#include "AbsBeamline/Degrader.h"
gsell's avatar
gsell committed
27 28 29 30
#include "AbsBeamline/Drift.h"
#include "AbsBeamline/SBend.h"
#include "AbsBeamline/RBend.h"
#include "AbsBeamline/Multipole.h"
31
#include "Structure/LossDataSink.h"
kraus's avatar
kraus committed
32
#include "Utilities/Options.h"
33
#include "Utilities/GeneralClassicException.h"
34 35
#include "Utilities/Util.h"
#include "Utilities/Timer.h"
adelmann's avatar
adelmann committed
36

37
#include "Utility/Inform.h"
38

snuverink_j's avatar
snuverink_j committed
39 40
#include <gsl/gsl_randist.h>

41
#include <algorithm>
snuverink_j's avatar
snuverink_j committed
42
#include <cmath>
gsell's avatar
gsell committed
43
#include <fstream>
44
#include <iostream>
45

46 47
#include <sys/time.h>

48 49
namespace {
    struct DegraderInsideTester: public InsideTester {
50
        explicit DegraderInsideTester(ElementBase* el) {
51 52 53
            deg_m = static_cast<Degrader*>(el);
        }
        virtual
54
        bool checkHit(const Vector_t& R) override {
55 56 57 58
            return deg_m->isInMaterial(R(2));
        }

    private:
59
        Degrader* deg_m;
60 61 62
    };

    struct CollimatorInsideTester: public InsideTester {
63
        explicit CollimatorInsideTester(ElementBase* el) {
64 65 66
            col_m = static_cast<CCollimator*>(el);
        }
        virtual
67
        bool checkHit(const Vector_t& R)  override {
68 69 70 71
            return col_m->checkPoint(R(0), R(1));
        }

    private:
72
        CCollimator* col_m;
73 74 75
    };

    struct FlexCollimatorInsideTester: public InsideTester {
76
        explicit FlexCollimatorInsideTester(ElementBase* el) {
77 78
            col_m = static_cast<FlexibleCollimator*>(el);
        }
79

80
        virtual
81
        bool checkHit(const Vector_t& R)  override {
82
            return col_m->isStopped(R);
83 84 85
        }

    private:
86
        FlexibleCollimator* col_m;
87 88 89
    };
}

90
ScatteringPhysics::ScatteringPhysics(const std::string& name,
91 92
                                     ElementBase* element,
                                     std::string& material,
ext-calvo_p's avatar
ext-calvo_p committed
93 94
                                     bool enableRutherford,
                                     double lowEnergyThr):
95
    ParticleMatterInteractionHandler(name, element),
96 97
    T_m(0.0),
    dT_m(0.0),
98 99
    mass_m(0.0),
    charge_m(0.0),
100
    material_m(material),
101
    hitTester_m(nullptr),
gsell's avatar
gsell committed
102 103
    Z_m(0),
    A_m(0.0),
104 105 106
    rho_m(0.0),
    X0_m(0.0),
    I_m(0.0),
107
    A1_c(0.0),
108 109 110 111
    A2_c(0.0),
    A3_c(0.0),
    A4_c(0.0),
    A5_c(0.0),
112 113 114 115 116
    B1_c(0.0),
    B2_c(0.0),
    B3_c(0.0),
    B4_c(0.0),
    B5_c(0.0),
117 118 119
    bunchToMatStat_m(0),
    stoppedPartStat_m(0),
    rediffusedStat_m(0),
120
    totalPartsInMat_m(0),
121 122
    Eavg_m(0.0),
    Emax_m(0.0),
123
    Emin_m(0.0),
ext-calvo_p's avatar
ext-calvo_p committed
124 125
    enableRutherford_m(enableRutherford),
    lowEnergyThr_m(lowEnergyThr)
126
{
adelmann's avatar
adelmann committed
127

adelmann's avatar
adelmann committed
128 129
    gsl_rng_env_setup();
    rGen_m = gsl_rng_alloc(gsl_rng_default);
130 131 132 133 134 135 136 137 138 139

    size_t mySeed = Options::seed;

    if (Options::seed == -1) {
        struct timeval tv;
        gettimeofday(&tv,0);
        mySeed = tv.tv_sec + tv.tv_usec;
    }

    gsl_rng_set(rGen_m, mySeed + Ippl::myNode());
adelmann's avatar
adelmann committed
140

141 142
    configureMaterialParameters();

143
    collshape_m = element_ref_m->getType();
144 145 146 147 148 149 150 151 152 153 154
    switch (collshape_m) {
    case ElementBase::DEGRADER:
        hitTester_m.reset(new DegraderInsideTester(element_ref_m));
        break;
    case ElementBase::CCOLLIMATOR:
        hitTester_m.reset(new CollimatorInsideTester(element_ref_m));
        break;
    case ElementBase::FLEXIBLECOLLIMATOR:
        hitTester_m.reset(new FlexCollimatorInsideTester(element_ref_m));
        break;
    default:
155
        throw GeneralClassicException("ScatteringPhysics::ScatteringPhysics",
156 157
                            "Unsupported element type");
    }
adelmann's avatar
adelmann committed
158

159
    lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(getName(), !Options::asciidump));
160

161 162
    DegraderApplyTimer_m   = IpplTimings::getTimer("DegraderApply");
    DegraderLoopTimer_m    = IpplTimings::getTimer("DegraderLoop");
163
    DegraderDestroyTimer_m = IpplTimings::getTimer("DegraderDestroy");
gsell's avatar
gsell committed
164 165
}

166
ScatteringPhysics::~ScatteringPhysics() {
167 168 169 170
    locParts_m.clear();
    lossDs_m->save();
    if (rGen_m)
        gsl_rng_free(rGen_m);
gsell's avatar
gsell committed
171 172
}

173

174 175
/// The material of the collimator
//  ------------------------------------------------------------------------
176
void  ScatteringPhysics::configureMaterialParameters() {
177

178 179 180 181 182 183
    auto material = Physics::Material::getMaterial(material_m);
    Z_m = material->getAtomicNumber();
    A_m = material->getAtomicMass();
    rho_m = material->getMassDensity();
    X0_m = material->getRadiationLength();
    I_m = material->getMeanExcitationEnergy();
184
    A1_c = material->getStoppingPowerFitCoefficients(Physics::Material::A1);
185 186 187 188
    A2_c = material->getStoppingPowerFitCoefficients(Physics::Material::A2);
    A3_c = material->getStoppingPowerFitCoefficients(Physics::Material::A3);
    A4_c = material->getStoppingPowerFitCoefficients(Physics::Material::A4);
    A5_c = material->getStoppingPowerFitCoefficients(Physics::Material::A5);
189 190 191 192 193
    B1_c = material->getStoppingPowerFitCoefficients(Physics::Material::B1);
    B2_c = material->getStoppingPowerFitCoefficients(Physics::Material::B2);
    B3_c = material->getStoppingPowerFitCoefficients(Physics::Material::B3);
    B4_c = material->getStoppingPowerFitCoefficients(Physics::Material::B4);
    B5_c = material->getStoppingPowerFitCoefficients(Physics::Material::B5);
adelmann's avatar
adelmann committed
194 195
}

196
void ScatteringPhysics::apply(PartBunchBase<double, 3>* bunch,
197
                              const std::pair<Vector_t, double>& boundingSphere) {
adelmann's avatar
adelmann committed
198
    IpplTimings::startTimer(DegraderApplyTimer_m);
199

200 201
    /*
      Particles that have entered material are flagged as Bin[i] == -1.
gsell's avatar
gsell committed
202

203
      Flagged particles are copied to a local structure within Scattering Physics locParts_m.
gsell's avatar
gsell committed
204

205
      Particles in that structure will be pushed in the material and either come
206
      back to the bunch or will be fully stopped in the material.
207
    */
208

kraus's avatar
kraus committed
209 210 211 212 213 214 215 216 217 218 219 220 221 222
    ParticleType pType = bunch->getPType();
    if (pType != ParticleType::PROTON   &&
        pType != ParticleType::DEUTERON &&
        pType != ParticleType::HMINUS   &&
        pType != ParticleType::MUON     &&
        pType != ParticleType::H2P      &&
        pType != ParticleType::ALPHA) {

        throw GeneralClassicException(
                "ScatteringPhysics::apply",
                "Particle " + bunch->getPTypeString() +
                " is not supported for scattering interactions!");
    }

adelmann's avatar
adelmann committed
223 224 225 226
    Eavg_m = 0.0;
    Emax_m = 0.0;
    Emin_m = 0.0;

227
    bunchToMatStat_m  = 0;
228
    rediffusedStat_m   = 0;
229
    stoppedPartStat_m = 0;
230
    totalPartsInMat_m   = 0;
231

frey_m's avatar
frey_m committed
232 233
    dT_m = bunch->getdT();
    T_m  = bunch->getT();
234 235
    mass_m   = bunch->getM();
    charge_m = bunch->getQ();
236

237
    bool onlyOneLoopOverParticles = ! (allParticleInMat_m);
238

239 240
    do {
        IpplTimings::startTimer(DegraderLoopTimer_m);
241 242
        push();
        setTimeStepForLeavingParticles();
243

244 245
        // if we are not looping copy newly arrived particles
        if (onlyOneLoopOverParticles) {
246
            copyFromBunch(bunch, boundingSphere);
247
        }
248
        addBackToBunch(bunch);
249

250
        computeInteraction(bunch);
251

252 253
        push();
        resetTimeStep();
254

255 256 257 258
        IpplTimings::stopTimer(DegraderLoopTimer_m);

        T_m += dT_m;              // update local time

259
        gatherStatistics();
260 261 262 263 264

        if (allParticleInMat_m) {
            onlyOneLoopOverParticles = (rediffusedStat_m > 0 || totalPartsInMat_m <= 0);
        } else {
            onlyOneLoopOverParticles = true;
265
        }
266
    } while (onlyOneLoopOverParticles == false);
snuverink_j's avatar
snuverink_j committed
267 268

    IpplTimings::stopTimer(DegraderApplyTimer_m);
269
}
270

271
void ScatteringPhysics::computeInteraction(PartBunchBase<double, 3>* bunch) {
272
    /*
273
        Do physics if
274
        -- correct type of particle
275
        -- particle not stopped (locParts_m[i].label != -1.0)
276

277 278
        Absorbed particle i: locParts_m[i].label = -1.0;
    */
kraus's avatar
kraus committed
279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301
    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;

            if (hitTester_m->checkHit(R)) {
                bool pdead = computeEnergyLoss(bunch, P, dt);
                if (!pdead) {
                    /*
                      Now scatter and transport particle in material.
                      The checkHit call just above will detect if the
                      particle is rediffused from the material into vacuum.
                    */

                    computeCoulombScattering(R, P, dt);
                } else {
                    // The particle is stopped in the material, set label to -1
                    locParts_m[i].label = -1.0;
                    ++ stoppedPartStat_m;
                    lossDs_m->addParticle(OpalParticle(locParts_m[i].IDincol,
                                                       R, P, T_m,
                                                       locParts_m[i].Qincol, locParts_m[i].Mincol));
302
                }
303
            }
304
        }
305
    }
kraus's avatar
kraus committed
306

307
    // delete absorbed particles
308 309
    deleteParticleFromLocalVector();
}
310

311
/// Energy Loss: using the Bethe-Bloch equation.
312
/// In low-energy region use Andersen-Ziegler fitting (only for protons and alpha)
313 314
/// Energy straggling: For relatively thick absorbers such that the number of collisions
/// is large, the energy loss distribution is shown to be Gaussian in form.
315 316 317
/// See Particle Physics Booklet, chapter 'Passage of particles through matter' or
/// Review of Particle Physics, DOI: 10.1103/PhysRevD.86.010001, page 329 ff
// -------------------------------------------------------------------------
kraus's avatar
kraus committed
318
bool ScatteringPhysics::computeEnergyLoss(PartBunchBase<double, 3>* bunch,
319
                                          Vector_t& P,
320 321 322
                                          const double deltat,
                                          bool includeFluctuations) const {

323
    ParticleType pType = bunch->getPType();
324
    constexpr double m2cm = 1e2;
325
    constexpr double eV2keV = 1e-3;
326
    constexpr double GeV2keV = 1e6;
327

328
    constexpr double massElectron_keV = Physics::m_e * GeV2keV;
329 330 331

    const double mass_keV = mass_m * eV2keV;

snuverink_j's avatar
snuverink_j committed
332
    constexpr double K = 4.0 * Physics::pi * Physics::Avo * Physics::r_e * m2cm * Physics::r_e * m2cm * massElectron_keV;
333 334 335 336

    double gamma = Util::getGamma(P);
    const double gammaSqr = std::pow(gamma, 2);
    const double betaSqr = 1.0 - 1.0 / gammaSqr;
snuverink_j's avatar
snuverink_j committed
337
    double beta = std::sqrt(betaSqr);
338
    double Ekin = (gamma - 1) * mass_keV;
339
    double dEdx = 0.0;
340
    double epsilon = 0.0;
341

342 343
    const double deltas = deltat * beta * Physics::c;
    const double deltasrho = deltas * m2cm * rho_m;
344

345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373
    const double massRatio = massElectron_keV / mass_keV;
    double Tmax = (2.0 * massElectron_keV * betaSqr * gammaSqr /
                  (std::pow(gamma + massRatio, 2) - (gammaSqr - 1.0)));

    if (pType != ParticleType::ALPHA) {

        if (Ekin >= 600 && Ekin < 1e7) {
            dEdx = (-K * std::pow(charge_m, 2) * Z_m / (A_m * betaSqr) *
                    (0.5 * std::log(2 * massElectron_keV * betaSqr * gammaSqr * Tmax / std::pow(I_m * eV2keV, 2)) - betaSqr));
        } else if (pType == ParticleType::PROTON && Ekin < 600) {
            constexpr double massProton_amu = Physics::m_p / Physics::amu;
            const double Ts = Ekin / massProton_amu;
            if (Ekin > 10) {
                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));
                epsilon = (epsilon_low * epsilon_high) / (epsilon_low + epsilon_high);
            } else if (Ekin > 1) {
                epsilon = A1_c * std::pow(Ts, 0.5);
            }
            dEdx = -epsilon / (1e18 * (A_m / Physics::Avo));
        } else {
            INFOMSG(level4 << "Particle energy out of the valid range "
                              "for energy loss calculation." << endl);
        }
    } else {
        if (Ekin > 10000 && Ekin < 1e6) {
            dEdx = (-K * std::pow(charge_m, 2) * Z_m / (A_m * betaSqr) *
                    (0.5 * std::log(2 * massElectron_keV * betaSqr * gammaSqr * Tmax / std::pow(I_m * eV2keV, 2)) - betaSqr));
        } else if (Ekin > 1 && Ekin <= 10000) {
kraus's avatar
kraus committed
374
            const double T = Ekin * 1e-3; //MeV
375 376 377 378 379 380 381 382
            const double epsilon_low = B1_c * std::pow(1000*T, B2_c);
            const double epsilon_high = (B3_c / T) * std::log(1 + (B4_c / T) + (B5_c * T));
            epsilon = (epsilon_low * epsilon_high) / (epsilon_low + epsilon_high);
            dEdx = -epsilon / (1e18 * (A_m / Physics::Avo));
        } else {
            INFOMSG(level4 << "Particle energy out of the valid range "
                              "for energy loss calculation." << endl);
        }
383
    }
384

385 386
    Ekin += deltasrho * dEdx;

387
    if (includeFluctuations) {
snuverink_j's avatar
snuverink_j committed
388
        double sigma_E = std::sqrt(K * massElectron_keV * rho_m * (Z_m / A_m) * deltas * m2cm);
389 390
        Ekin += gsl_ran_gaussian(rGen_m, sigma_E);
    }
391

392
    gamma = Ekin / mass_keV + 1.0;
snuverink_j's avatar
snuverink_j committed
393
    beta = std::sqrt(1.0 - 1.0 / std::pow(gamma, 2));
394
    P = gamma * beta * P / euclidean_norm(P);
adelmann's avatar
adelmann committed
395

ext-calvo_p's avatar
ext-calvo_p committed
396
    bool stopped = (Ekin < lowEnergyThr_m * 1e3 || dEdx > 0);
397 398
    return stopped;
}
adelmann's avatar
adelmann committed
399

400

401
// Implement the rotation in 2 dimensions here
402 403 404
// For details see: J. Beringer et al. (Particle Data Group),
// "Passage of particles through matter", Phys. Rev. D 86, 010001 (2012)
// -------------------------------------------------------------------------
405
void  ScatteringPhysics::applyRotation(Vector_t& P,
406
                                       Vector_t& R,
407 408 409 410 411
                                       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);

412 413
    R(0) = R(0) + shift * std::cos(Psixz);
    R(2) = R(2) - shift * std::sin(Psixz);
414 415 416

    // Apply the rotation about the random angle thetacou
    double Px = P(0);
417 418
    P(0) =  Px * std::cos(thetacou) + P(2) * std::sin(thetacou);
    P(2) = -Px * std::sin(thetacou) + P(2) * std::cos(thetacou);
gsell's avatar
gsell committed
419 420
}

421
void ScatteringPhysics::applyRandomRotation(Vector_t& P, double theta0) {
gsell's avatar
gsell committed
422

snuverink_j's avatar
snuverink_j committed
423
    double thetaru = 2.5 / std::sqrt(gsl_rng_uniform(rGen_m)) * 2.0 * theta0;
424
    double phiru = Physics::two_pi * gsl_rng_uniform(rGen_m);
gsell's avatar
gsell committed
425

snuverink_j's avatar
snuverink_j committed
426
    double normPtrans = std::sqrt(P(0) * P(0) + P(1) * P(1));
427
    double Theta = std::atan(normPtrans / std::abs(P(2)));
428 429
    double CosT = std::cos(Theta);
    double SinT = std::sin(Theta);
430

431 432 433
    Vector_t X(std::cos(phiru)*std::sin(thetaru),
               std::sin(phiru)*std::sin(thetaru),
               std::cos(thetaru));
434 435 436 437
    X *= euclidean_norm(P);

    Vector_t W(-P(1), P(0), 0.0);
    W = W / normPtrans;
gsell's avatar
gsell committed
438

439 440
    // Rodrigues' formula for a rotation about W by angle Theta
    P = X * CosT + cross(W, X) * SinT + W * dot(W, X) * (1.0 - CosT);
441 442
}

443 444
/// Coulomb Scattering: Including Multiple Coulomb Scattering and large angle Rutherford Scattering.
/// Using the distribution given in Classical Electrodynamics, by J. D. Jackson.
gsell's avatar
gsell committed
445
//--------------------------------------------------------------------------
446
void  ScatteringPhysics::computeCoulombScattering(Vector_t& R,
447
                                                  Vector_t& P,
448 449
                                                  double dt) {

snuverink_j's avatar
snuverink_j committed
450
    constexpr double sqrtThreeInv = 0.57735026918962576451; // sqrt(1.0 / 3.0)
451 452
    const double normP = euclidean_norm(P);
    const double gammaSqr = std::pow(normP, 2) + 1.0;
snuverink_j's avatar
snuverink_j committed
453
    const double beta = std::sqrt(1.0 - 1.0 / gammaSqr);
454
    const double deltas = dt * beta * Physics::c;
455 456 457
    const double theta0 = (13.6e6 / (beta * normP * mass_m) *
                           charge_m * std::sqrt(deltas / X0_m) *
                           (1.0 + 0.038 * std::log(deltas / X0_m)));
458 459 460 461 462 463 464 465 466 467 468 469 470 471

    double phi = Physics::two_pi * gsl_rng_uniform(rGen_m);
    for (unsigned int i = 0; i < 2; ++ i) {
        CoordinateSystemTrafo randomTrafo(R, Quaternion(cos(phi), 0, 0, sin(phi)));
        P = randomTrafo.rotateTo(P);
        R = Vector_t(0.0); // corresponds to randomTrafo.transformTo(R);

        double z1 = gsl_ran_gaussian(rGen_m, 1.0);
        double z2 = gsl_ran_gaussian(rGen_m, 1.0);

        while(std::abs(z2) > 3.5) {
            z1 = gsl_ran_gaussian(rGen_m, 1.0);
            z2 = gsl_ran_gaussian(rGen_m, 1.0);
        }
472

473 474 475
        double thetacou = z2 * theta0;
        double shift = (z1 * sqrtThreeInv + z2) * deltas * theta0 * 0.5;
        applyRotation(P, R, shift, thetacou);
gsell's avatar
gsell committed
476

477 478
        P = randomTrafo.rotateFrom(P);
        R = randomTrafo.transformFrom(R);
479

480
        phi += 0.5 * Physics::pi;
adelmann's avatar
adelmann committed
481
    }
482 483 484 485 486

    if (enableRutherford_m &&
        gsl_rng_uniform(rGen_m) < 0.0047) {

        applyRandomRotation(P, theta0);
gsell's avatar
gsell committed
487 488 489
    }
}

490
void ScatteringPhysics::addBackToBunch(PartBunchBase<double, 3>* bunch) {
adelmann's avatar
adelmann committed
491

492 493
    const size_t nL = locParts_m.size();
    if (nL == 0) return;
adelmann's avatar
adelmann committed
494

495 496
    unsigned int numLocalParticles = bunch->getLocalNum();
    const double elementLength = element_ref_m->getElementLength();
497

498
    for (size_t i = 0; i < nL; ++ i) {
499
        Vector_t& R = locParts_m[i].Rincol;
500 501 502 503 504 505 506 507 508 509 510 511 512

        if (R[2] >= elementLength) {

            bunch->createWithID(locParts_m[i].IDincol);

            /*
              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;
            bunch->Q[numLocalParticles]   = locParts_m[i].Qincol;
kraus's avatar
kraus committed
513
            bunch->M[numLocalParticles]   = locParts_m[i].Mincol;
514 515 516 517 518 519 520 521 522 523 524 525 526 527 528
            bunch->Bf[numLocalParticles]  = 0.0;
            bunch->Ef[numLocalParticles]  = 0.0;
            bunch->dt[numLocalParticles]  = dT_m;

            /*
              This particle is back to the bunch, by set
              ting the label to -1.0
              the particle will be deleted.
            */
            locParts_m[i].label = -1.0;

            ++ rediffusedStat_m;
            ++ numLocalParticles;
        }
    }
529

530
    // delete particles that went to the bunch
531
    deleteParticleFromLocalVector();
532 533
}

534
void ScatteringPhysics::copyFromBunch(PartBunchBase<double, 3>* bunch,
535
                                      const std::pair<Vector_t, double>& boundingSphere)
536
{
frey_m's avatar
frey_m committed
537
    const size_t nL = bunch->getLocalNum();
538 539
    if (nL == 0) return;

540
    IpplTimings::startTimer(DegraderDestroyTimer_m);
541 542 543
    double zmax = boundingSphere.first(2) + boundingSphere.second;
    double zmin = boundingSphere.first(2) - boundingSphere.second;
    if (zmax < 0.0 || zmin > element_ref_m->getElementLength()) {
544 545 546 547
        IpplTimings::stopTimer(DegraderDestroyTimer_m);
        return;
    }

adelmann's avatar
adelmann committed
548
    size_t ne = 0;
549
    std::set<size_t> partsToDel;
550
    for (size_t i = 0; i < nL; ++ i) {
frey_m's avatar
frey_m committed
551
        if ((bunch->Bin[i] == -1 || bunch->Bin[i] == 1) &&
552
            hitTester_m->checkHit(bunch->R[i]))
553
        {
554
            // adjust the time step for those particles that enter the material
snuverink_j's avatar
snuverink_j committed
555
            // such that it corresponds to the time needed to reach the current
556 557 558
            // location form the edge of the material. Only use this time step
            // for the computation of the interaction with the material, not for
            // the integration of the particles. This will ensure that the momenta
559
            // of all particles are reduced by approximately the same amount in
560 561 562 563 564 565 566 567
            // computeEnergyLoss.
            double betaz = bunch->P[i](2) / Util::getGamma(bunch->P[i]);
            double stepWidth = betaz * Physics::c * bunch->dt[i];
            double tau = bunch->R[i](2) / stepWidth;

            PAssert_LE(tau, 1.0);
            PAssert_GE(tau, 0.0);

568 569
            PART x;
            x.localID      = i;
570
            x.DTincol      = bunch->dt[i] * tau;
frey_m's avatar
frey_m committed
571 572 573 574 575
            x.IDincol      = bunch->ID[i];
            x.Binincol     = bunch->Bin[i];
            x.Rincol       = bunch->R[i];
            x.Pincol       = bunch->P[i];
            x.Qincol       = bunch->Q[i];
kraus's avatar
kraus committed
576
            x.Mincol       = bunch->M[i];
frey_m's avatar
frey_m committed
577 578
            x.Bfincol      = bunch->Bf[i];
            x.Efincol      = bunch->Ef[i];
snuverink_j's avatar
snuverink_j committed
579
            x.label        = 0;            // alive in matter
580 581 582 583 584

            locParts_m.push_back(x);
            ne++;
            bunchToMatStat_m++;

585
            partsToDel.insert(i);
586
        }
587
    }
588

589
    for (auto it = partsToDel.begin(); it != partsToDel.end(); ++ it) {
590
        bunch->destroy(1, *it);
591 592
    }

593
    if (ne > 0) {
frey_m's avatar
frey_m committed
594
        bunch->performDestroy(true);
595
    }
596

597
    IpplTimings::stopTimer(DegraderDestroyTimer_m);
598 599
}

600
void ScatteringPhysics::print(Inform &msg) {
601
    Inform::FmtFlags_t ff = msg.flags();
kraus's avatar
kraus committed
602

603
    if (totalPartsInMat_m > 0 ||
604 605
        bunchToMatStat_m  > 0 ||
        rediffusedStat_m  > 0 ||
606 607
        stoppedPartStat_m > 0) {

608
        OPALTimer::Timer time;
609
        msg << level2
610
            << "--- ScatteringPhysics - Name: " << name_m << "\n"
611
            << "Material: " << material_m << " - Element: " << element_ref_m->getName() << "\n"
612 613 614 615
            << "Particle Statistics @ " << time.time() << "\n"
            << std::setw(21) << "entered: " << Util::toStringWithThousandSep(bunchToMatStat_m) << "\n"
            << std::setw(21) << "rediffused: " << Util::toStringWithThousandSep(rediffusedStat_m) << "\n"
            << std::setw(21) << "stopped: " << Util::toStringWithThousandSep(stoppedPartStat_m) << "\n"
616
            << std::setw(21) << "total in material: " << Util::toStringWithThousandSep(totalPartsInMat_m)
617
            << endl;
618
    }
619 620

    msg.flags(ff);
621
}
gsell's avatar
gsell committed
622

623
bool ScatteringPhysics::stillActive() {
624
    return totalPartsInMat_m != 0;
625
}
626

627 628 629 630
namespace {
    bool myCompF(PART x, PART y) {
      return x.label > y.label;
    }
631 632
}

633
void ScatteringPhysics::deleteParticleFromLocalVector() {
634
    /*
kraus's avatar
kraus committed
635
      the particle to be deleted (label < 0) are all
636
      at the end of the vector.
637
    */
638
    sort(locParts_m.begin(), locParts_m.end(), myCompF);
639

adelmann's avatar
adelmann committed
640 641
    // find start of particles to delete
    std::vector<PART>::iterator inv = locParts_m.begin();
642

643
    for (; inv != locParts_m.end(); ++inv) {
adelmann's avatar
adelmann committed
644 645 646
        if ((*inv).label == -1)
            break;
    }
647 648
    locParts_m.erase(inv, locParts_m.end());
    locParts_m.resize(inv - locParts_m.begin());
adelmann's avatar
adelmann committed
649 650 651 652 653 654

    // update statistics
    if (locParts_m.size() > 0) {
        Eavg_m /= locParts_m.size();
        Emin_m /= locParts_m.size();
        Emax_m /= locParts_m.size();
655
    }
gsell's avatar
gsell committed
656
}
657

658
void ScatteringPhysics::push() {
659
    for (size_t i = 0; i < locParts_m.size(); ++i) {
660 661
        Vector_t& R  = locParts_m[i].Rincol;
        Vector_t& P  = locParts_m[i].Pincol;
662 663 664 665 666 667 668 669 670 671 672 673
        double gamma = Util::getGamma(P);

        R += 0.5 * dT_m * Physics::c * P / gamma;
    }
}

// adjust the time step for those particles that will rediffuse to
// vacuum such that it corresponds to the time needed to reach the
// end of the material. Only use this time step for the computation
// of the interaction with the material, not for the integration of
// the particles. This will ensure that the momenta of all particles
// are reduced by  approcimately the same amount in computeEnergyLoss.
674
void ScatteringPhysics::setTimeStepForLeavingParticles() {
675 676 677
    const double elementLength = element_ref_m->getElementLength();

    for (size_t i = 0; i < locParts_m.size(); ++i) {
678 679 680
        Vector_t& R  = locParts_m[i].Rincol;
        Vector_t& P  = locParts_m[i].Pincol;
        double& dt   = locParts_m[i].DTincol;
681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697
        double gamma = Util::getGamma(P);
        Vector_t stepLength = dT_m * Physics::c * P / gamma;

        if (R(2) < elementLength &&
            R(2) + stepLength(2) > elementLength) {
            // particle is likely to leave material
            double distance = elementLength - R(2);
            double tau = distance / stepLength(2);

            PAssert_LE(tau, 1.0);
            PAssert_GE(tau, 0.0);

            dt *= tau;
        }
    }
}

698
void ScatteringPhysics::resetTimeStep() {
699
    for (size_t i = 0; i < locParts_m.size(); ++i) {
700
        double& dt = locParts_m[i].DTincol;
701 702 703 704 705
        dt = dT_m;
    }

}

706

707
void ScatteringPhysics::gatherStatistics() {
708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723

    unsigned int locPartsInMat;
    locPartsInMat = locParts_m.size();

    constexpr unsigned short numItems = 4;
    unsigned int partStatistics[numItems] = {locPartsInMat,
                                             bunchToMatStat_m,
                                             rediffusedStat_m,
                                             stoppedPartStat_m};

    allreduce(&(partStatistics[0]), numItems, std::plus<unsigned int>());

    totalPartsInMat_m = partStatistics[0];
    bunchToMatStat_m = partStatistics[1];
    rediffusedStat_m = partStatistics[2];
    stoppedPartStat_m = partStatistics[3];
724
}