TrackRun.cpp

Go to the documentation of this file.

//
// Class TrackRun
//   The RUN command.
//
// Copyright (c) 200x - 2023, 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/>.
//
#include "Track/TrackRun.h"

#include "AbstractObjects/BeamSequence.h"
#include "AbstractObjects/ObjectFunction.h"
#include "AbstractObjects/OpalData.h"

#include "Algorithms/ParallelCyclotronTracker.h"
#include "Algorithms/ParallelTTracker.h"
#include "Algorithms/ThickTracker.h"
#include "Algorithms/Tracker.h"

#include "Attributes/Attributes.h"

#include "Beamlines/TBeamline.h"

#include "BasicActions/Option.h"

#include "Distribution/Distribution.h"

#include "Physics/Physics.h"
#include "Physics/Units.h"

#include "Track/Track.h"

#include "Utilities/OpalException.h"

#include "Structure/Beam.h"
#include "Structure/BoundaryGeometry.h"
#include "Structure/DataSink.h"
#include "Structure/FieldSolver.h"
#include "Structure/H5PartWrapper.h"
#include "Structure/H5PartWrapperForPC.h"
#include "Structure/H5PartWrapperForPT.h"

#include "OPALconfig.h"
#include "changes.h"

#include <boost/assign.hpp>

#include <cmath>
#include <fstream>
#include <iomanip>


extern Inform *gmsg;

std::shared_ptr<Tracker> TrackRun::itsTracker_m = nullptr;

namespace {
    // The attributes of class TrackRun.
    enum {
        METHOD,           // Tracking method to use.
        TURNS,            // The number of turns to be tracked.
        MBMODE,           // The working way for multi-bunch mode for OPAL-cycl: "FORCE" or "AUTO"
        PARAMB,           // The control parameter for "AUTO" mode of multi-bunch,
        MB_ETA,           // The scale parameter for binning in multi-bunch mode
        MB_BINNING,       // The binning type in multi-bunch mode
        BEAM,             // The beam to track
        FIELDSOLVER,      // The field solver attached
        BOUNDARYGEOMETRY, // The boundary geometry
        DISTRIBUTION,     // The particle distribution
        TRACKBACK,
        SIZE
    };
}

const std::string TrackRun::defaultDistribution("DISTRIBUTION");

const boost::bimap<TrackRun::RunMethod, std::string> TrackRun::stringMethod_s =
    boost::assign::list_of<const boost::bimap<TrackRun::RunMethod, std::string>::relation>
    (RunMethod::PARALLELT,  "PARALLEL-T")
    (RunMethod::CYCLOTRONT, "CYCLOTRON-T")
    (RunMethod::THICK,      "THICK");


TrackRun::TrackRun():
    Action(SIZE, "RUN",
           "The \"RUN\" sub-command tracks the defined particles through "
           "the given lattice."),
    dist_m(nullptr),
    fieldSolver_m(nullptr),
    dataSink_m(nullptr),
    phaseSpaceSink_m(nullptr),
    isFollowupTrack_m(false),
    method_m(RunMethod::NONE),
    macromass_m(0.0),
    macrocharge_m(0.0) {
    itsAttr[METHOD] = Attributes::makePredefinedString
                      ("METHOD", "Name of tracking algorithm to use.",
                       {"THICK", "OPAL-T", "PARALLEL-T", "OPAL-CYCL", "CYCLOTRON-T"});
    
    itsAttr[TURNS] = Attributes::makeReal
        ("TURNS", "Number of turns to be tracked; Number of neighboring bunches to be tracked in cyclotron.", 1.0);

    itsAttr[MBMODE] = Attributes::makePredefinedString
        ("MBMODE", "The working way for multi-bunch mode for OPAL-cycl.",
         {"FORCE", "AUTO"}, "FORCE");

    itsAttr[PARAMB] = Attributes::makeReal
        ("PARAMB", "Control parameter to define when to start multi-bunch mode, only available in \"AUTO\" mode.", 5.0);

    itsAttr[MB_ETA] = Attributes::makeReal
        ("MB_ETA", "The scale parameter for binning in multi-bunch mode.", 0.01);

    itsAttr[MB_BINNING] = Attributes::makePredefinedString
        ("MB_BINNING", "Type of energy binning in multi-bunch mode.",
         {"GAMMA_BINNING", "BUNCH_BINNING"}, "GAMMA_BINNING");

    itsAttr[BEAM] = Attributes::makeString
        ("BEAM", "Name of beam.");

    itsAttr[FIELDSOLVER] = Attributes::makeString
        ("FIELDSOLVER", "Field solver to be used.");

    itsAttr[BOUNDARYGEOMETRY] = Attributes::makeString
        ("BOUNDARYGEOMETRY", "Boundary geometry to be used NONE (default).", "NONE");

    itsAttr[DISTRIBUTION] = Attributes::makeStringArray
        ("DISTRIBUTION", "List of particle distributions to be used.");

    itsAttr[TRACKBACK] = Attributes::makeBool
        ("TRACKBACK", "Track in reverse direction, default: false.", false);

    registerOwnership(AttributeHandler::SUB_COMMAND);
    opalData_m = OpalData::getInstance();
}


TrackRun::TrackRun(const std::string& name, TrackRun* parent):
    Action(name, parent),
    dist_m(nullptr),
    fieldSolver_m(nullptr),
    dataSink_m(nullptr),
    phaseSpaceSink_m(nullptr),
    isFollowupTrack_m(false),
    method_m(RunMethod::NONE),
    macromass_m(0.0),
    macrocharge_m(0.0) {
    opalData_m = OpalData::getInstance();
}


TrackRun::~TrackRun() {
    delete phaseSpaceSink_m;
}


TrackRun* TrackRun::clone(const std::string& name) {
    return new TrackRun(name, this);
}


void TrackRun::execute() {
    const int currentVersion = ((OPAL_VERSION_MAJOR * 100) + OPAL_VERSION_MINOR) * 100;
    if (Options::version < currentVersion) {
        unsigned int fileVersion = Options::version / 100;
        bool newerChanges = false;
        for (auto it = Versions::changes.begin(); it != Versions::changes.end(); ++ it) {
            if (it->first > fileVersion) {
                newerChanges = true;
                break;
            }
        }
        if (newerChanges) {
            Inform errorMsg("Error");
            errorMsg << "\n******************** V E R S I O N   M I S M A T C H ***********************\n" << endl;
            for (auto it = Versions::changes.begin(); it != Versions::changes.end(); ++ it) {
                if (it->first > fileVersion) {
                    errorMsg << it->second << endl;
                }
            }
            errorMsg << "\n"
                     << "* Make sure you do understand these changes and adjust your input file \n"
                     << "* accordingly. Then add\n"
                     << "* OPTION, VERSION = " << currentVersion << ";\n"
                     << "* to your input file. " << endl;
            errorMsg << "\n****************************************************************************\n" << endl;
            throw OpalException("TrackRun::execute", "Version mismatch");
        }
    }

    isFollowupTrack_m = opalData_m->hasBunchAllocated();
    if (!itsAttr[DISTRIBUTION] && !isFollowupTrack_m) {
        throw OpalException("TrackRun::execute",
                            "\"DISTRIBUTION\" must be set in \"RUN\" command.");
    }
    if (!itsAttr[FIELDSOLVER]) {
        throw OpalException("TrackRun::execute",
                            "\"FIELDSOLVER\" must be set in \"RUN\" command.");
    }
    if (!itsAttr[BEAM]) {
        throw OpalException("TrackRun::execute",
                            "\"BEAM\" must be set in \"RUN\" command.");
    }

    // Get algorithm to use.
    setRunMethod();
    switch (method_m) {
       case RunMethod::THICK: {
            setupThickTracker();
            break;
        }
        case RunMethod::PARALLELT: {
            setupTTracker();
            break;
        }
        case RunMethod::CYCLOTRONT: {
            setupCyclotronTracker();
            break;
        }
        default: {
            throw OpalException("TrackRun::execute",
                                "Unknown \"METHOD\" for the \"RUN\" command");
        }
    }

    if (method_m == RunMethod::THICK) {
        int turns = int(std::round(Attributes::getReal(itsAttr[TURNS])));

        // Track for the all but last turn.
        for (int turn = 1; turn < turns; ++turn) {
            itsTracker_m->execute();
        }
        // Track the last turn.
        itsTracker_m->execute();

    } else {
        itsTracker_m->execute();

        opalData_m->setRestartRun(false);
    }

    opalData_m->bunchIsAllocated();
}

void TrackRun::setRunMethod() {
    if (!itsAttr[METHOD]) {
        throw OpalException("TrackRun::setRunMethod",
                            "The attribute \"METHOD\" isn't set for the \"RUN\" command");
    } else {
        auto it = stringMethod_s.right.find(Attributes::getString(itsAttr[METHOD]));
        if (it != stringMethod_s.right.end()) {
            method_m = it->second;
        }
    }
}

std::string TrackRun::getRunMethodName() const {
    return stringMethod_s.left.at(method_m);
}

void TrackRun::setupThickTracker() {
    if (isFollowupTrack_m) {
        Track::block->bunch->setLocalTrackStep(0);
    }

    Beam* beam = Beam::find(Attributes::getString(itsAttr[BEAM]));

    setBoundaryGeometry();

    setupFieldsolver();

    initPhaseSpaceSink();

    macrocharge_m = setDistributionParallelT(beam);

    *gmsg << *this << endl;

    Track::block->bunch->setdT(Track::block->dT.front());
    Track::block->bunch->dtScInit_m = Track::block->dtScInit;
    Track::block->bunch->deltaTau_m = Track::block->deltaTau;

    if (!isFollowupTrack_m && !opalData_m->inRestartRun()) {
        Track::block->bunch->setT(Track::block->t0_m);
    }

    if (Track::block->bunch->getIfBeamEmitting()) {
        Track::block->bunch->setChargeZeroPart(macrocharge_m);
    } else {
        Track::block->bunch->setCharge(macrocharge_m);
    }

    // set coupling constant
    double coefE = 1.0 / (4 * Physics::pi * Physics::epsilon_0);
    Track::block->bunch->setCouplingConstant(coefE);

    // statistical data are calculated (rms, eps etc.)
    Track::block->bunch->calcBeamParameters();

    initDataSink();

    if (!isFollowupTrack_m) {
        *gmsg << *dist_m << endl;
    }

    if (Track::block->bunch->getTotalNum() > 0) {
        double spos = /*Track::block->bunch->get_sPos() +*/ Track::block->zstart;
        auto &zstop = Track::block->zstop;
        auto &timeStep = Track::block->localTimeSteps;
        auto &dT = Track::block->dT;

        unsigned int i = 0;
        while (i + 1 < zstop.size() && zstop[i + 1] < spos) {
            ++ i;
        }

        zstop.erase(zstop.begin(), zstop.begin() + i);
        timeStep.erase(timeStep.begin(), timeStep.begin() + i);
        dT.erase(dT.begin(), dT.begin() + i);

        Track::block->bunch->setdT(dT.front());
    } else {
        Track::block->zstart = 0.0;
    }

    *gmsg << *beam << endl;
    *gmsg << *fieldSolver_m   << endl;

    itsTracker_m.reset(new ThickTracker(*Track::block->use->fetchLine(),
                                  Track::block->bunch, *beam, *dataSink_m, Track::block->reference,
                                  false, false, Track::block->localTimeSteps,
                                  Track::block->zstart, Track::block->zstop, Track::block->dT,
                                  Track::block->truncOrder));
}


void TrackRun::setupTTracker(){
    opalData_m->setInOPALTMode();

    if (isFollowupTrack_m) {
        Track::block->bunch->setLocalTrackStep(0);
    }

    Beam* beam = Beam::find(Attributes::getString(itsAttr[BEAM]));
    Track::block->bunch->setBeamFrequency(beam->getFrequency() * Units::MHz2Hz);
    Track::block->bunch->setPType(beam->getParticleName());

    setBoundaryGeometry();

    setupFieldsolver();

    initPhaseSpaceSink();

    macrocharge_m = setDistributionParallelT(beam);
    macromass_m   = beam->getMassPerParticle();

    *gmsg << *this << endl;

    Track::block->bunch->setdT(Track::block->dT.front());
    Track::block->bunch->dtScInit_m = Track::block->dtScInit;
    Track::block->bunch->deltaTau_m = Track::block->deltaTau;

    if (!isFollowupTrack_m && !opalData_m->inRestartRun()) {
        Track::block->bunch->setT(Track::block->t0_m);
    }

    if (Track::block->bunch->getIfBeamEmitting()) {
        Track::block->bunch->setChargeZeroPart(macrocharge_m);
        Track::block->bunch->setMassZeroPart(macromass_m);
    } else {
        Track::block->bunch->setCharge(macrocharge_m);
        Track::block->bunch->setMass(macromass_m);
    }
    // set coupling constant
    double coefE = 1.0 / (4 * Physics::pi * Physics::epsilon_0);
    Track::block->bunch->setCouplingConstant(coefE);

    // statistical data are calculated (rms, eps etc.)
    Track::block->bunch->calcBeamParameters();

    initDataSink();

    if (!isFollowupTrack_m) {
        *gmsg << std::scientific;
        *gmsg << *dist_m << endl;
    }

    if (Track::block->bunch->getTotalNum() > 0) {
        double spos = Track::block->zstart;
        auto& zstop = Track::block->zstop;
        auto it = Track::block->dT.begin();

        unsigned int i = 0;
        while (i + 1 < zstop.size() && zstop[i + 1] < spos) {
            ++ i;
            ++ it;
        }

        Track::block->bunch->setdT(*it);
    } else {
        Track::block->zstart = 0.0;
    }

    *gmsg << *beam << endl;
    *gmsg << *fieldSolver_m   << endl;

    // findPhasesForMaxEnergy();

    itsTracker_m.reset(new ParallelTTracker(*Track::block->use->fetchLine(),
                                      Track::block->bunch,
                                      *dataSink_m,
                                      Track::block->reference,
                                      false,
                                      Attributes::getBool(itsAttr[TRACKBACK]),
                                      Track::block->localTimeSteps,
                                      Track::block->zstart,
                                      Track::block->zstop,
                                      Track::block->dT));
}

void TrackRun::setupCyclotronTracker(){
    opalData_m->setInOPALCyclMode();

    Beam* beam = Beam::find(Attributes::getString(itsAttr[BEAM]));

    setBoundaryGeometry();

    setupFieldsolver();

    Track::block->bunch->setPType(beam->getParticleName());
    Track::block->bunch->POrigin = ParticleOrigin::REGULAR;

    std::vector<std::string> distr_str = Attributes::getStringArray(itsAttr[DISTRIBUTION]);
    if (distr_str.size() == 0) {
        dist_m = Distribution::find(defaultDistribution);
    } else {
        dist_m = Distribution::find(distr_str.at(0));
    }

    // multi-bunch parameters
    const int specifiedNumBunch = int(std::abs(std::round(Attributes::getReal(itsAttr[TURNS]))));
    const double mbPara         = Attributes::getReal(itsAttr[PARAMB]);
    const std::string mbMode    = Attributes::getString(itsAttr[MBMODE]);
    const double mbEta          = Attributes::getReal(itsAttr[MB_ETA]);
    const std::string mbBinning = Attributes::getString(itsAttr[MB_BINNING]);

    initPhaseSpaceSink();

    if (beam->getNumberOfParticles() < 3 || beam->getCurrent() == 0.0) {
        macrocharge_m = beam->getCharge() * Physics::q_e;
        macromass_m   = beam->getMass();
        Track::block->bunch->setDistribution(dist_m,
                                             beam->getNumberOfParticles(),
                                             beam->getCurrent(),
                                             *Track::block->use->fetchLine());

    } else {
        macrocharge_m = beam->getChargePerParticle();
        macromass_m   = beam->getMassPerParticle();

        if (!isFollowupTrack_m) {
            if (!opalData_m->inRestartRun()) {
                Track::block->bunch->setDistribution(dist_m,
                                                     beam->getNumberOfParticles(),
                                                     beam->getCurrent(),
                                                     *Track::block->use->fetchLine());

            } else {
                dist_m->doRestartOpalCycl(Track::block->bunch,
                                        beam->getNumberOfParticles(),
                                        opalData_m->getRestartStep(),
                                        specifiedNumBunch,
                                        phaseSpaceSink_m);
            }
        }
    }
    Track::block->bunch->setMass(macromass_m); // set the Mass per macro-particle, [GeV/c^2]
    Track::block->bunch->setCharge(macrocharge_m);  // set the charge per macro-particle, [C]

    Track::block->bunch->setdT(1.0 / (Track::block->stepsPerTurn * beam->getFrequency() * Units::MHz2Hz));
    Track::block->bunch->setStepsPerTurn(Track::block->stepsPerTurn);

    // set coupling constant
    double coefE = 1.0 / (4 * Physics::pi * Physics::epsilon_0);
    Track::block->bunch->setCouplingConstant(coefE);

    // statistical data are calculated (rms, eps etc.)
    Track::block->bunch->calcBeamParameters();

    initDataSink(specifiedNumBunch);

    itsTracker_m.reset(new ParallelCyclotronTracker(*Track::block->use->fetchLine(),
                                              Track::block->bunch, *dataSink_m, Track::block->reference,
                                              false, false, Track::block->localTimeSteps.front(),
                                              Track::block->timeIntegrator,
                                              specifiedNumBunch, mbEta, mbPara, mbMode, mbBinning));

    ParallelCyclotronTracker* cyclTracker = dynamic_cast<ParallelCyclotronTracker*>(itsTracker_m.get());

    if (opalData_m->inRestartRun()) {
        H5PartWrapperForPC *h5pw = static_cast<H5PartWrapperForPC*>(phaseSpaceSink_m);
        cyclTracker->setBeGa(h5pw->getMeanMomentum());

        cyclTracker->setPr(h5pw->getReferencePr());
        cyclTracker->setPt(h5pw->getReferencePt());
        cyclTracker->setPz(h5pw->getReferencePz());

        cyclTracker->setR(h5pw->getReferenceR());
        cyclTracker->setTheta(h5pw->getReferenceT());
        cyclTracker->setZ(h5pw->getReferenceZ());

        // The following is for restarts in local frame
        cyclTracker->setPhi(h5pw->getAzimuth());
        cyclTracker->setPsi(h5pw->getElevation());
        cyclTracker->setPreviousH5Local(h5pw->getPreviousH5Local());

        if ( specifiedNumBunch > 1 ) {
            cyclTracker->setLastDumpedStep(opalData_m->getRestartStep());
        }
    }

    // statistical data are calculated (rms, eps etc.)
    Track::block->bunch->calcBeamParameters();

    *gmsg << *this << endl;
    *gmsg << *dist_m << endl;
    *gmsg << *beam << endl;
    *gmsg << *fieldSolver_m   << endl;
}

void TrackRun::setupFieldsolver() {
    fieldSolver_m = FieldSolver::find(Attributes::getString(itsAttr[FIELDSOLVER]));

    if (fieldSolver_m->getFieldSolverType() != FieldSolverType::NONE) {
        size_t numGridPoints = fieldSolver_m->getMX()*fieldSolver_m->getMY()*fieldSolver_m->getMT(); // total number of gridpoints
        Beam* beam = Beam::find(Attributes::getString(itsAttr[BEAM]));
        size_t numParticles = beam->getNumberOfParticles();

        if (!opalData_m->inRestartRun() && numParticles < numGridPoints
            && fieldSolver_m->getFieldSolverType() != FieldSolverType::SAAMG // in SPIRAL/SAAMG we're meshing the whole domain -DW
            && fieldSolver_m->getFieldSolverType() != FieldSolverType::P3M //In P3M with one-one mapping grid points can be less than particles
            && !Options::amr)
        {
            throw OpalException("TrackRun::setupFieldsolver()",
                                "The number of simulation particles (" + std::to_string(numParticles) + ") \n" +
                                "is smaller than the number of gridpoints (" + std::to_string(numGridPoints) + ").\n" +
                                "Please increase the number of particles or reduce the size of the mesh.\n");
        }

        opalData_m->addProblemCharacteristicValue("MX", fieldSolver_m->getMX());
        opalData_m->addProblemCharacteristicValue("MY", fieldSolver_m->getMY());
        opalData_m->addProblemCharacteristicValue("MT", fieldSolver_m->getMT());
    }

    fieldSolver_m->initCartesianFields();
    Track::block->bunch->setSolver(fieldSolver_m);
    if (fieldSolver_m->hasPeriodicZ()) {
        Track::block->bunch->setBCForDCBeam();
    } else {
        Track::block->bunch->setBCAllOpen();
    }
}

void TrackRun::initPhaseSpaceSink() {
    const std::string h5FileName = opalData_m->getInputBasename() + std::string(".h5");
    if (opalData_m->isInOPALCyclMode()) {
        if (opalData_m->inRestartRun()) {
            phaseSpaceSink_m = new H5PartWrapperForPC(h5FileName,
                                                      opalData_m->getRestartStep(),
                                                      opalData_m->getRestartFileName(),
                                                      H5_O_WRONLY);
        } else if (isFollowupTrack_m) {
            phaseSpaceSink_m = new H5PartWrapperForPC(h5FileName,
                                                      -1,
                                                      h5FileName,
                                                      H5_O_WRONLY);
        } else if (Options::enableHDF5) {
            phaseSpaceSink_m = new H5PartWrapperForPC(h5FileName,
                                                      H5_O_WRONLY);
        }
    } else {
        if (opalData_m->inRestartRun()) {
            phaseSpaceSink_m = new H5PartWrapperForPT(h5FileName,
                                                      opalData_m->getRestartStep(),
                                                      opalData_m->getRestartFileName(),
                                                      H5_O_WRONLY);
        } else if (isFollowupTrack_m) {
            phaseSpaceSink_m = new H5PartWrapperForPT(h5FileName,
                                                      -1,
                                                      h5FileName,
                                                      H5_O_WRONLY);
        } else if (Options::enableHDF5) {
            phaseSpaceSink_m = new H5PartWrapperForPT(h5FileName,
                                                      H5_O_WRONLY);
        }
    }
}


void TrackRun::initDataSink(const int& numBunch) {
    if (!opalData_m->inRestartRun()) {
        if (!opalData_m->hasDataSinkAllocated()) {
            opalData_m->setDataSink(new DataSink(phaseSpaceSink_m, false, numBunch));
        } else {
            dataSink_m = opalData_m->getDataSink();
            dataSink_m->changeH5Wrapper(phaseSpaceSink_m);
        }
    } else {
        opalData_m->setDataSink(new DataSink(phaseSpaceSink_m, true, numBunch));
    }
    dataSink_m = opalData_m->getDataSink();
}

void TrackRun::setBoundaryGeometry() {
    if (Attributes::getString(itsAttr[BOUNDARYGEOMETRY]) != "NONE") {
        // Ask the dictionary if BoundaryGeometry is allocated.
        // If it is allocated use the allocated BoundaryGeometry
        if (!opalData_m->hasGlobalGeometry()) {
            const std::string geomDescriptor = Attributes::getString(itsAttr[BOUNDARYGEOMETRY]);
            BoundaryGeometry* bg = BoundaryGeometry::find(geomDescriptor)->clone(geomDescriptor);
            opalData_m->setGlobalGeometry(bg);
        }
    }
}


double TrackRun::setDistributionParallelT(Beam* beam) {
    /*
     * Distribution(s) can be set via a single distribution or a list
     * (array) of distributions. If an array is defined the first in the
     * list is treated as the primary distribution. All others are added to
     * it to create the full distribution.
     */
    std::vector<std::string> distributionArray
        = Attributes::getStringArray(itsAttr[DISTRIBUTION]);
    const size_t numberOfDistributions = distributionArray.size();

    if (numberOfDistributions == 0) {
        dist_m = Distribution::find(defaultDistribution);
    } else {
        dist_m = Distribution::find(distributionArray.at(0));
        dist_m->setNumberOfDistributions(numberOfDistributions);

        if (numberOfDistributions > 1) {
            *gmsg << endl
                  << "---------------------------------" << endl
                  << "Found more than one distribution:" << endl << endl;
            *gmsg << "Main Distribution" << endl
                  << "---------------------------------" << endl
                  << distributionArray.at(0) << endl << endl
                  << "Secondary Distribution(s)" << endl
                  << "---------------------------------" << endl;

            for (size_t i = 1; i < numberOfDistributions; ++ i) {
                Distribution *distribution = Distribution::find(distributionArray.at(i));
                distribution->setNumberOfDistributions(numberOfDistributions);
                distrs_m.push_back(distribution);

                *gmsg << distributionArray.at(i) << endl;
            }
            *gmsg << endl
                  << "---------------------------------" << endl << endl;
        }
    }

    /*
     * Initialize distributions.
     */
    size_t numberOfParticles = beam->getNumberOfParticles();
    if (!isFollowupTrack_m) {
        if (!opalData_m->inRestartRun()) {
            /*
             * Here we are not doing a restart run
             * and we do not have a bunch already allocated.
             */
            Track::block->bunch->setDistribution(dist_m,
                                                 distrs_m,
                                                 numberOfParticles);

            /*
             * If this is an injected beam (rather than an emitted beam), we
             * make sure it doesn't have any particles at z < 0.
             */

            opalData_m->setGlobalPhaseShift(0.5 * dist_m->getTEmission() + dist_m->getEmissionTimeShift());
        } else {
            /*
             * Read in beam from restart file.
             */
            dist_m->doRestartOpalT(Track::block->bunch, numberOfParticles, opalData_m->getRestartStep(), phaseSpaceSink_m);
        }
    }

    // Return charge per macroparticle.
    return beam->getChargePerParticle();
}

Inform& TrackRun::print(Inform& os) const {
    os << endl;
    os << "* ************* T R A C K  R U N *************************************************** " << endl;
    if (!isFollowupTrack_m) {
        os << "* Selected Tracking Method == " << getRunMethodName() << ", NEW TRACK" << '\n'
           << "* ********************************************************************************** " << '\n';
    } else {
        os << "* Selected Tracking Method == " << getRunMethodName() << ", FOLLOWUP TRACK" << '\n'
           << "* ********************************************************************************** " << '\n';
    }
    os << "* Phase space dump frequency    = " << Options::psDumpFreq << '\n'
       << "* Statistics dump frequency     = " << Options::statDumpFreq << " w.r.t. the time step." << '\n'
       << "* DT                            = " << Track::block->dT.front() << " [s]\n"
       << "* MAXSTEPS                      = " << Track::block->localTimeSteps.front() << '\n'
       << "* Mass of simulation particle   = " << macromass_m << " [GeV/c^2]" << '\n'
       << "* Charge of simulation particle = " << macrocharge_m << " [C]" << '\n';
    if (method_m == RunMethod::CYCLOTRONT) {
        os << "* Number of neighbour bunches   = " << int(std::abs(std::round(Attributes::getReal(itsAttr[TURNS])))) << '\n'
           << "* STEPSPERTURN                  = " << Track::block->stepsPerTurn << '\n';
    }
    os << "* ********************************************************************************** ";
    return os;
}

std::shared_ptr<Tracker> TrackRun::getTracker() {
    return itsTracker_m;
}