ParallelSliceTracker.cpp 23 KB
Newer Older
kraus's avatar
kraus committed
1
#include "Algorithms/ParallelSliceTracker.h"
2 3 4
#include "Algorithms/OrbitThreader.h"
#include "Algorithms/IndexMap.h"
#include "Algorithms/CavityAutophaser.h"
gsell's avatar
gsell committed
5

6
#include "AbstractObjects/OpalData.h"
gsell's avatar
gsell committed
7
#include "Beamlines/Beamline.h"
8
#include "Beamlines/FlaggedBeamline.h"
gsell's avatar
gsell committed
9
#include "Distribution/Distribution.h"
10
#include "Lines/Sequence.h"
gsell's avatar
gsell committed
11
#include "Utilities/Timer.h"
12
#include "Utilities/Util.h"
13
#include "AbsBeamline/SpecificElementVisitor.h"
gsell's avatar
gsell committed
14

15 16 17
#include <memory>
#include <string>

gsell's avatar
gsell committed
18 19 20
class PartData;

ParallelSliceTracker::ParallelSliceTracker(const Beamline &beamline,
21
        const PartData &reference, bool revBeam, bool revTrack):
gsell's avatar
gsell committed
22
    Tracker(beamline, reference, revBeam, revTrack) {
23 24

    CoordinateSystemTrafo labToRef(beamline.getOrigin3D(),
25
                                   beamline.getInitialDirection());
26 27
    referenceToLabCSTrafo_m = labToRef.inverted();

gsell's avatar
gsell committed
28 29 30 31
}


ParallelSliceTracker::ParallelSliceTracker(const Beamline &beamline,
32 33 34 35 36 37 38 39 40
                                           EnvelopeBunch &bunch,
                                           DataSink &ds,
                                           const PartData &reference,
                                           bool revBeam,
                                           bool revTrack,
                                           const std::vector<unsigned long long> &maxSteps,
                                           double zstart,
                                           const std::vector<double> &zstop,
                                           const std::vector<double> &dt):
gsell's avatar
gsell committed
41
    Tracker(beamline, reference, revBeam, revTrack),
42 43
    itsOpalBeamline_m(),
    zstart_m(zstart) {
gsell's avatar
gsell committed
44

45 46
    itsBunch_m      = &bunch;
    itsDataSink_m   = &ds;
gsell's avatar
gsell committed
47

48
    CoordinateSystemTrafo labToRef(beamline.getOrigin3D(),
49
                                   beamline.getInitialDirection());
50 51 52 53 54 55 56 57 58 59 60
    referenceToLabCSTrafo_m = labToRef.inverted();

    for (std::vector<unsigned long long>::const_iterator it = maxSteps.begin(); it != maxSteps.end(); ++ it) {
        localTrackSteps_m.push(*it);
    }
    for (std::vector<double>::const_iterator it = dt.begin(); it != dt.end(); ++ it) {
        dtAllTracks_m.push(*it);
    }
    for (std::vector<double>::const_iterator it = zstop.begin(); it != zstop.end(); ++ it) {
        zStop_m.push(*it);
    }
gsell's avatar
gsell committed
61 62 63

    timeIntegrationTimer1_m  = IpplTimings::getTimer("Time integration1");
    timeIntegrationTimer2_m  = IpplTimings::getTimer("Time integration2");
64 65 66
    timeFieldEvaluation_m    = IpplTimings::getTimer("Field evaluation");
    BinRepartTimer_m         = IpplTimings::getTimer("Time of Binary repart.");
    WakeFieldTimer_m         = IpplTimings::getTimer("Time of Wake Field calc.");
67

68 69
    cavities_m = itsOpalBeamline_m.getElementByType(ElementBase::RFCAVITY);
    auto travelingwaves = itsOpalBeamline_m.getElementByType(ElementBase::TRAVELINGWAVE);
70
    cavities_m.merge(travelingwaves, ClassicField::SortAsc);
gsell's avatar
gsell committed
71 72 73
}


74
ParallelSliceTracker::~ParallelSliceTracker()
75 76 77 78 79 80
{ }


void ParallelSliceTracker::visitBeamline(const Beamline &bl) { // borrowed from ParallelTTracker
    const FlaggedBeamline* fbl = static_cast<const FlaggedBeamline*>(&bl);
    if (fbl->getRelativeFlag()) {
81
        OpalBeamline stash(fbl->getOrigin3D(), fbl->getInitialDirection());
82 83 84 85 86 87 88 89 90 91
        stash.swap(itsOpalBeamline_m);
        fbl->iterate(*this, false);
        itsOpalBeamline_m.prepareSections();
        itsOpalBeamline_m.compute3DLattice();
        stash.merge(itsOpalBeamline_m);
        stash.swap(itsOpalBeamline_m);
    } else {
        fbl->iterate(*this, false);
    }
}
gsell's avatar
gsell committed
92 93


94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110
void ParallelSliceTracker::saveCavityPhases() { // borrowed from ParallelTTracker
    itsDataSink_m->storeCavityInformation();
}

void ParallelSliceTracker::restoreCavityPhases() { // borrowed from ParallelTTracker
    typedef std::vector<MaxPhasesT>::iterator iterator_t;

    if (OpalData::getInstance()->hasPriorTrack() ||
        OpalData::getInstance()->inRestartRun()) {
        iterator_t it = OpalData::getInstance()->getFirstMaxPhases();
        iterator_t end = OpalData::getInstance()->getLastMaxPhases();
        for (; it < end; ++ it) {
            updateRFElement((*it).first, (*it).second);
        }
    }
}

111 112 113 114 115
/*
 * The maximum phase is added to the nominal phase of
 * the element. This is done on all nodes except node 0 where
 * the Autophase took place.
 */
116 117
void ParallelSliceTracker::updateRFElement(std::string elName, double maxPhase) {
    double phase = 0.0;
gsell's avatar
gsell committed
118

119
    for (FieldList::iterator fit = cavities_m.begin(); fit != cavities_m.end(); ++fit) {
120
        if ((*fit).getElement()->getName() == elName) {
121
            if ((*fit).getElement()->getType() == ElementBase::TRAVELINGWAVE) {
122
                phase  =  static_cast<TravelingWave *>((*fit).getElement().get())->getPhasem();
123
                static_cast<TravelingWave *>((*fit).getElement().get())->setPhasem(phase + maxPhase);
gsell's avatar
gsell committed
124
            } else {
125
                phase  = static_cast<RFCavity *>((*fit).getElement().get())->getPhasem();
126
                static_cast<RFCavity *>((*fit).getElement().get())->setPhasem(phase + maxPhase);
gsell's avatar
gsell committed
127
            }
128 129

            break;
gsell's avatar
gsell committed
130 131 132 133
        }
    }
}

134 135 136 137 138

/*
 * All RF-Elements gets updated, where the phiShift is the
 * global phase shift in units of seconds.
 */
139 140
void ParallelSliceTracker::printRFPhases() {
    Inform msg("ParallelSliceTracker ");
gsell's avatar
gsell committed
141

142
    FieldList &cl = cavities_m;
gsell's avatar
gsell committed
143

144 145
    const double RADDEG = 180.0 / Physics::pi;
    const double globalTimeShift = OpalData::getInstance()->getGlobalPhaseShift();
146 147

    msg << "\n-------------------------------------------------------------------------------------\n";
gsell's avatar
gsell committed
148

149
    for (FieldList::iterator it = cl.begin(); it != cl.end(); ++it) {
150 151 152 153 154
        std::shared_ptr<Component> element(it->getElement());
        std::string name = element->getName();
        double frequency;
        double phase;

155
        if (element->getType() == ElementBase::TRAVELINGWAVE) {
156
            phase = static_cast<TravelingWave *>(element.get())->getPhasem();
kraus's avatar
kraus committed
157
            frequency = static_cast<TravelingWave *>(element.get())->getFrequencym();
gsell's avatar
gsell committed
158
        } else {
159
            phase = static_cast<RFCavity *>(element.get())->getPhasem();
kraus's avatar
kraus committed
160
            frequency = static_cast<RFCavity *>(element.get())->getFrequencym();
gsell's avatar
gsell committed
161
        }
162 163 164 165 166

        msg << (it == cl.begin()? "": "\n")
            << name
            << ": phi = phi_nom + phi_maxE + global phase shift = " << (phase - globalTimeShift * frequency) * RADDEG << " degree, "
            << "(global phase shift = " << -globalTimeShift *frequency *RADDEG << " degree) \n";
gsell's avatar
gsell committed
167
    }
168

169
    msg << "-------------------------------------------------------------------------------------\n"
kraus's avatar
kraus committed
170
        << endl;
gsell's avatar
gsell committed
171 172 173 174
}

void ParallelSliceTracker::execute() {

175 176 177
    Inform msg("ParallelSliceTracker", *gmsg);
    BorisPusher pusher(itsReference);
    const double globalTimeShift = itsBunch_m->doEmission()? OpalData::getInstance()->getGlobalPhaseShift(): 0.0;
gsell's avatar
gsell committed
178

179
    currentSimulationTime_m = itsBunch_m->getT();
gsell's avatar
gsell committed
180

181 182 183 184
    setTime();
    setLastStep();

    prepareSections();
185
    handleAutoPhasing();
gsell's avatar
gsell committed
186

187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231
    std::queue<double> timeStepSizes(dtAllTracks_m);
    std::queue<unsigned long long> numSteps(localTrackSteps_m);
    double minTimeStep = timeStepSizes.front();
    unsigned long long totalNumSteps = 0;
    while (timeStepSizes.size() > 0) {
        if (minTimeStep > timeStepSizes.front()) {
            totalNumSteps = std::ceil(totalNumSteps * minTimeStep / timeStepSizes.front());
            minTimeStep = timeStepSizes.front();
        }
        totalNumSteps += std::ceil(numSteps.front() * timeStepSizes.front() / minTimeStep);

        numSteps.pop();
        timeStepSizes.pop();
    }

    itsOpalBeamline_m.activateElements();

    if ( OpalData::getInstance()->hasPriorTrack() ||
         OpalData::getInstance()->inRestartRun()) {

        referenceToLabCSTrafo_m = itsBunch_m->toLabTrafo_m;
        RefPartR_m = referenceToLabCSTrafo_m.transformFrom(itsBunch_m->RefPartR_m);
        RefPartP_m = referenceToLabCSTrafo_m.rotateFrom(itsBunch_m->RefPartP_m);

        pathLength_m = itsBunch_m->get_sPos();
        zstart_m = pathLength_m;

        restoreCavityPhases();
    } else {
        RefPartR_m = Vector_t(0.0);
        RefPartP_m = euclidean_norm(itsBunch_m->get_pmean_Distribution()) * Vector_t(0, 0, 1);

        if (itsBunch_m->getTotalNum() > 0) {
            if (!itsOpalBeamline_m.containsSource()) {
                RefPartP_m = itsReference.getP() / itsBunch_m->getM() * Vector_t(0, 0, 1);
            }

            if (zstart_m > pathLength_m) {
                findStartPosition(pusher);
            }

            itsBunch_m->set_sPos(pathLength_m);
        }
    }

gsell's avatar
gsell committed
232
    Vector_t rmin, rmax;
233 234 235 236 237 238 239 240 241 242 243 244 245 246 247
    itsBunch_m->get_bounds(rmin, rmax);
    OrbitThreader oth(itsReference,
                      referenceToLabCSTrafo_m.transformTo(RefPartR_m),
                      referenceToLabCSTrafo_m.rotateTo(RefPartP_m),
                      pathLength_m,
                      0.0,//-rmin(2),
                      itsBunch_m->getT(),
                      minTimeStep,
                      totalNumSteps,
                      zStop_m.back() + 2 * rmax(2),
                      itsOpalBeamline_m);

    oth.execute();

    saveCavityPhases();
gsell's avatar
gsell committed
248

249
    setTime();
gsell's avatar
gsell committed
250

251 252 253 254 255 256 257
    double t = itsBunch_m->getT() - globalTimeShift;
    itsBunch_m->setT(t);

    itsBunch_m->RefPartR_m = referenceToLabCSTrafo_m.transformTo(RefPartR_m);
    itsBunch_m->RefPartP_m = referenceToLabCSTrafo_m.rotateTo(RefPartP_m);

    *gmsg << level1 << *itsBunch_m << endl;
gsell's avatar
gsell committed
258

259 260 261 262 263
    unsigned long long step = itsBunch_m->getGlobalTrackStep();
    OPALTimer::Timer myt1;
    *gmsg << "Track start at: " << myt1.time() << ", t= " << Util::getTimeString(t) << "; "
          << "zstart at: " << Util::getLengthString(pathLength_m)
          << endl;
gsell's avatar
gsell committed
264

265
    // prepareEmission();
266

267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283
    *gmsg << level1
          << "Executing ParallelSliceTracker, initial dt= " << Util::getTimeString(itsBunch_m->getdT()) << ";\n"
          << "max integration steps " << getMaxSteps(localTrackSteps_m) << ", next step= " << step << endl
          << "the mass is: " << itsReference.getM() * 1e-6 << " MeV, its charge: "
          << itsReference.getQ() << " [e]" << endl;


    // setOptionalVariables();
    *gmsg << std::scientific;
    itsBunch_m->toLabTrafo_m = referenceToLabCSTrafo_m;
    while (localTrackSteps_m.size() > 0) {
        localTrackSteps_m.front() += step;
        dtCurrentTrack_m = dtAllTracks_m.front();
        selectDT();

        for (; step < localTrackSteps_m.front(); ++step) {

284
            globalEOL_m = false;
285 286 287 288 289 290

            switchElements();
            computeExternalFields(oth);

            //reduce(&globalEOL_m, &globalEOL_m, OpBitwiseOrAssign());
            //reduce(&globalEOL_m, &globalEOL_m + 1, &globalEOL_m, OpBitwiseAndAssign());
291
            allreduce(&globalEOL_m, 1, std::logical_and<bool>());
292 293 294 295 296 297 298 299 300 301 302

            computeSpaceChargeFields();
            timeIntegration();

            if (step % 1000 + 1 == 1000) {
                msg << level1;
            } else if (step % 100 + 1 == 100) {
                msg << level2;
            } else {
                msg << level3;
            }
303 304 305 306 307

            msg << " Step " << step << " at " << itsBunch_m->zAvg() << " [m] t= "
                << itsBunch_m->getT() << " [s] E=" << itsBunch_m->Eavg() * 1e-6
                << " [MeV]" << endl;

308 309
            currentSimulationTime_m += itsBunch_m->getdT();
            itsBunch_m->setT(currentSimulationTime_m);
310

311 312
            dumpStats(step);

313 314 315 316 317 318 319 320 321
            if (hasEndOfLineReached())
                break;

            double gamma = itsBunch_m->Eavg() / (Physics::m_e * 1e9) + 1.0;
            double beta = sqrt(1.0 - 1.0 / (gamma * gamma));

            double driftPerTimeStep = itsBunch_m->getdT() * Physics::c * beta;
            if (std::abs(zStop_m.front() - itsBunch_m->zAvg()) < 0.5 * driftPerTimeStep)
                localTrackSteps_m.front() = step;
322
        }
323

324
        if (globalEOL_m)
325 326 327 328 329
            break;

        dtAllTracks_m.pop();
        localTrackSteps_m.pop();
        zStop_m.pop();
gsell's avatar
gsell committed
330 331
    }

332 333
    OpalData::getInstance()->setLastStep(step);
    writeLastStepPhaseSpace(step, itsBunch_m->get_sPos());
334
    itsOpalBeamline_m.switchElementsOff();
335 336 337 338 339 340 341 342 343 344 345 346
    msg << "done executing ParallelSliceTracker" << endl;
}


void ParallelSliceTracker::timeIntegration() {

    IpplTimings::startTimer(timeIntegrationTimer1_m);
    itsBunch_m->timeStep(itsBunch_m->getdT());
    IpplTimings::stopTimer(timeIntegrationTimer1_m);
}


347
void ParallelSliceTracker::handleAutoPhasing() {
348

349
    if (Options::autoPhase == 0) return;
gsell's avatar
gsell committed
350

351
    if(!OpalData::getInstance()->inRestartRun()) {
352
        itsDataSink_m->storeCavityInformation();
gsell's avatar
gsell committed
353 354
    }

355 356 357 358 359 360
    auto it = OpalData::getInstance()->getFirstMaxPhases();
    auto end = OpalData::getInstance()->getLastMaxPhases();
    for (; it != end; ++ it) {
        updateRFElement((*it).first, (*it).second);
    }
    printRFPhases();
361
}
gsell's avatar
gsell committed
362 363


364
void ParallelSliceTracker::prepareSections() {
gsell's avatar
gsell committed
365

366
    itsBeamline_m.accept(*this);
367
    itsOpalBeamline_m.prepareSections();
gsell's avatar
gsell committed
368

369 370 371
    itsOpalBeamline_m.compute3DLattice();
    itsOpalBeamline_m.save3DLattice();
    itsOpalBeamline_m.save3DInput();
372
}
gsell's avatar
gsell committed
373 374


375
void ParallelSliceTracker::computeExternalFields(OrbitThreader &oth) {
gsell's avatar
gsell committed
376

377
    IpplTimings::startTimer(timeFieldEvaluation_m);
gsell's avatar
gsell committed
378

379
    Vector_t externalE, externalEThis, externalB, externalBThis, KR, KRThis, KT, KTThis;
gsell's avatar
gsell committed
380

381 382 383 384 385 386 387 388 389 390 391 392
    const unsigned int localNum = itsBunch_m->getLocalNum();
    // bool locPartOutOfBounds = false, globPartOutOfBounds = false;
    Vector_t rmin, rmax;
    itsBunch_m->get_bounds(rmin, rmax);
    IndexMap::value_t elements;

    try {
        elements = oth.query(pathLength_m + 0.5 * (rmax(2) + rmin(2)), rmax(2) - rmin(2));
    } catch(IndexMap::OutOfBounds &e) {
        globalEOL_m = true;
        return;
    }
gsell's avatar
gsell committed
393

394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430
    const double time = itsBunch_m->getT() + 0.5 * itsBunch_m->getdT();
    const IndexMap::value_t::const_iterator end = elements.end();

    for (unsigned int i = 0; i < localNum; ++ i) {

        externalB     = Vector_t(0.0);
        externalE     = Vector_t(0.0);
        KR            = Vector_t(0.0);
        KT            = Vector_t(0.0);

        IndexMap::value_t::const_iterator it = elements.begin();
        for (; it != end; ++ it) {
            const CoordinateSystemTrafo toLocal = (itsOpalBeamline_m.getMisalignment((*it)) *
                                                   itsOpalBeamline_m.getCSTrafoLab2Local((*it)));
            const CoordinateSystemTrafo fromLocal = toLocal.inverted();

            //FIXME: why not x=y=0.0?
            // for (unsigned int k = 0; k < 100; ++ k) {
            externalEThis = Vector_t(0.0);
            externalBThis = Vector_t(0.0);
            KRThis        = Vector_t(0.0);
            KTThis        = Vector_t(0.0);

            Vector_t pos = itsBunch_m->getR(i);
            Vector_t mom = itsBunch_m->getP(i);

            itsBunch_m->setR(i, toLocal.transformTo(pos));
            itsBunch_m->setP(i, toLocal.rotateTo(mom));

            if ((*it)->apply(itsBunch_m->getR(i),
                             itsBunch_m->getP(i),
                             time,
                             externalEThis,
                             externalBThis)) {
                itsBunch_m->Bin[i] = -1;
                continue;
            }
gsell's avatar
gsell committed
431

432 433
            (*it)->addKR(i, time, KRThis);
            (*it)->addKT(i, time, KTThis);
gsell's avatar
gsell committed
434

435 436
            externalE += fromLocal.rotateTo(externalEThis);
            externalB += fromLocal.rotateTo(externalBThis);
gsell's avatar
gsell committed
437

438 439
            KR += fromLocal.rotateTo(KRThis);
            KT += fromLocal.rotateTo(KTThis);
gsell's avatar
gsell committed
440

441 442 443
            itsBunch_m->setR(i, pos);
            itsBunch_m->setP(i, mom);
        }
gsell's avatar
gsell committed
444

445 446
        itsBunch_m->setExternalFields(i, externalE, externalB, KR, KT);
    }
gsell's avatar
gsell committed
447

448 449 450 451 452 453 454 455
    IpplTimings::stopTimer(timeFieldEvaluation_m);
}


void ParallelSliceTracker::computeSpaceChargeFields() {

    itsBunch_m->computeSpaceCharge();
}
gsell's avatar
gsell committed
456 457


458
void ParallelSliceTracker::dumpStats(long long step) {
gsell's avatar
gsell committed
459

snuverink_j's avatar
snuverink_j committed
460 461 462
    double sposRef = itsBunch_m->get_sPos();
    if (step != 0 && (step % Options::psDumpFreq == 0 || step % Options::statDumpFreq == 0))
        writePhaseSpace(step, sposRef);
463
}
gsell's avatar
gsell committed
464

465

466
void ParallelSliceTracker::switchElements(double /*scaleMargin*/) {
467

snuverink_j's avatar
snuverink_j committed
468
    double margin = 1.0;
469

snuverink_j's avatar
snuverink_j committed
470
    currentSimulationTime_m = itsBunch_m->getT();
471
    itsOpalBeamline_m.switchElements(itsBunch_m->zTail() - margin,
snuverink_j's avatar
snuverink_j committed
472 473
                                      itsBunch_m->zHead() + margin,
                                      itsBunch_m->Eavg() * 1e-6);
474 475 476 477 478
}


void ParallelSliceTracker::setLastStep() {

479
    // unsigned long long step = 0;
480

481
    // if (OpalData::getInstance()->inRestartRun()) {
482

483 484 485 486
    //     int prevDumpFreq = OpalData::getInstance()->getRestartDumpFreq();
    //     step = OpalData::getInstance()->getRestartStep() * prevDumpFreq + 1;
    //     maxSteps_m += step;
    // } else {
487

488 489 490
    //     step = OpalData::getInstance()->getLastStep() + 1;
    //     maxSteps_m += step;
    // }
gsell's avatar
gsell committed
491

492
    // OpalData::getInstance()->setLastStep(step);
493
}
gsell's avatar
gsell committed
494 495


496
void ParallelSliceTracker::setTime() {
gsell's avatar
gsell committed
497

498 499 500 501 502 503 504 505
    if (OpalData::getInstance()->inRestartRun())
        currentSimulationTime_m = itsBunch_m->getT();
    else
        currentSimulationTime_m = itsBunch_m->getT();
}


bool ParallelSliceTracker::hasEndOfLineReached() {
506 507
    reduce(&globalEOL_m, &globalEOL_m + 1, &globalEOL_m, OpBitwiseAndAssign());
    return globalEOL_m;
508
}
509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666

void ParallelSliceTracker::findStartPosition(const BorisPusher &pusher) { // borrowed from ParallelTTracker

    double t = 0.0;
    itsBunch_m->setT(t);

    dtCurrentTrack_m = dtAllTracks_m.front();
    selectDT();

    if (Util::getEnergy(RefPartP_m, itsBunch_m->getM()) < 1e-3) {
        double gamma = 0.1 / itsBunch_m->getM() + 1.0;
        RefPartP_m = sqrt(std::pow(gamma, 2) - 1) * Vector_t(0, 0, 1);
    }

    while (true) {
        autophaseCavities(pusher);

        t += itsBunch_m->getdT();
        itsBunch_m->setT(t);

        Vector_t oldR = RefPartR_m;
        updateReferenceParticle(pusher);
        pathLength_m += euclidean_norm(RefPartR_m - oldR);

        if (pathLength_m > zStop_m.front()) {
            if (localTrackSteps_m.size() == 0) return;

            dtAllTracks_m.pop();
            localTrackSteps_m.pop();
            zStop_m.pop();

            selectDT();
        }

        double speed = euclidean_norm(RefPartP_m) * Physics::c / sqrt(dot(RefPartP_m, RefPartP_m) + 1);
        if (std::abs(pathLength_m - zstart_m) <=  0.5 * itsBunch_m->getdT() * speed) {
            double tau = (pathLength_m - zstart_m) / speed;

            t += tau;
            itsBunch_m->setT(t);

            RefPartR_m /= (Physics::c * tau);
            pusher.push(RefPartR_m, RefPartP_m, tau);
            RefPartR_m *= (Physics::c * tau);

            pathLength_m = zstart_m;

            CoordinateSystemTrafo update(RefPartR_m,
                                         getQuaternion(RefPartP_m, Vector_t(0, 0, 1)));
            referenceToLabCSTrafo_m = referenceToLabCSTrafo_m * update.inverted();

            RefPartR_m = update.transformTo(RefPartR_m);
            RefPartP_m = update.rotateTo(RefPartP_m);

            return;
        }
    }
}

void ParallelSliceTracker::updateReferenceParticle(const BorisPusher &pusher) { // borrowed from ParallelTTracker
    //static size_t step = 0;
    const double dt = std::min(itsBunch_m->getT(), itsBunch_m->getdT());
    const double scaleFactor = Physics::c * dt;
    Vector_t Ef(0.0), Bf(0.0);
    // Vector_t oldR = RefPartR_m;

    RefPartR_m /= scaleFactor;
    pusher.push(RefPartR_m, RefPartP_m, dt);
    RefPartR_m *= scaleFactor;

    IndexMap::value_t elements = itsOpalBeamline_m.getElements(referenceToLabCSTrafo_m.transformTo(RefPartR_m));
    IndexMap::value_t::const_iterator it = elements.begin();
    const IndexMap::value_t::const_iterator end = elements.end();

    for (; it != end; ++ it) {
        CoordinateSystemTrafo refToLocalCSTrafo = itsOpalBeamline_m.getCSTrafoLab2Local((*it)) * referenceToLabCSTrafo_m;

        Vector_t localR = refToLocalCSTrafo.transformTo(RefPartR_m);
        Vector_t localP = refToLocalCSTrafo.rotateTo(RefPartP_m);
        Vector_t localE(0.0), localB(0.0);

        if ((*it)->applyToReferenceParticle(localR,
                                            localP,
                                            itsBunch_m->getT() - 0.5 * dt,
                                            localE,
                                            localB)) {
            *gmsg << level1 << "The reference particle hit an element" << endl;
            globalEOL_m = true;
        }

        Ef += refToLocalCSTrafo.rotateFrom(localE);
        Bf += refToLocalCSTrafo.rotateFrom(localB);
    }

    pusher.kick(RefPartR_m, RefPartP_m, Ef, Bf, dt);

    RefPartR_m /= scaleFactor;
    pusher.push(RefPartR_m, RefPartP_m, dt);
    RefPartR_m *= scaleFactor;
    //++ step;
}

void ParallelSliceTracker::autophaseCavities(const BorisPusher &pusher) { // borrowed from ParallelTTracker

    double t = itsBunch_m->getT();
    Vector_t nextR = RefPartR_m / (Physics::c * itsBunch_m->getdT());
    pusher.push(nextR, RefPartP_m, itsBunch_m->getdT());
    nextR *= Physics::c * itsBunch_m->getdT();

    auto elementSet = itsOpalBeamline_m.getElements(referenceToLabCSTrafo_m.transformTo(nextR));
    for (auto element: elementSet) {
        if (element->getType() == ElementBase::TRAVELINGWAVE) {
            const TravelingWave *TWelement = static_cast<const TravelingWave *>(element.get());
            if (!TWelement->getAutophaseVeto()) {
                RefPartR_m = referenceToLabCSTrafo_m.transformTo(RefPartR_m);
                RefPartP_m = referenceToLabCSTrafo_m.rotateTo(RefPartP_m);
                CavityAutophaser ap(itsReference, element);
                ap.getPhaseAtMaxEnergy(itsOpalBeamline_m.transformToLocalCS(element, RefPartR_m),
                                       itsOpalBeamline_m.rotateToLocalCS(element, RefPartP_m),
                                       t, itsBunch_m->getdT());
                RefPartR_m = referenceToLabCSTrafo_m.transformFrom(RefPartR_m);
                RefPartP_m = referenceToLabCSTrafo_m.rotateFrom(RefPartP_m);
            }

        } else if (element->getType() == ElementBase::RFCAVITY) {
            const RFCavity *RFelement = static_cast<const RFCavity *>(element.get());
            if (!RFelement->getAutophaseVeto()) {
                RefPartR_m = referenceToLabCSTrafo_m.transformTo(RefPartR_m);
                RefPartP_m = referenceToLabCSTrafo_m.rotateTo(RefPartP_m);
                CavityAutophaser ap(itsReference, element);
                ap.getPhaseAtMaxEnergy(itsOpalBeamline_m.transformToLocalCS(element, RefPartR_m),
                                       itsOpalBeamline_m.rotateToLocalCS(element, RefPartP_m),
                                       t, itsBunch_m->getdT());
                RefPartR_m = referenceToLabCSTrafo_m.transformFrom(RefPartR_m);
                RefPartP_m = referenceToLabCSTrafo_m.rotateFrom(RefPartP_m);
            }
        }
    }
}

void ParallelSliceTracker::selectDT() { // borrowed from ParallelTTracker

    if (itsBunch_m->getIfBeamEmitting()) {
        double dt = itsBunch_m->getEmissionDeltaT();
        itsBunch_m->setdT(dt);
    } else {
        double dt = dtCurrentTrack_m;
        itsBunch_m->setdT(dt);
    }
}

void ParallelSliceTracker::changeDT() { // borrowed from ParallelTTracker
    selectDT();
    const unsigned int localNum = itsBunch_m->getLocalNum();
    for (unsigned int i = 0; i < localNum; ++ i) {
        itsBunch_m->dt[i] = itsBunch_m->getdT();
    }
}