ParallelSliceTracker.cpp 23.3 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 157
            phase = static_cast<TravelingWave *>(element.get())->getPhasem();
	    frequency = static_cast<TravelingWave *>(element.get())->getFrequencym();
gsell's avatar
gsell committed
158
        } else {
159 160
            phase = static_cast<RFCavity *>(element.get())->getPhasem();
	    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 170
    msg << "-------------------------------------------------------------------------------------\n"
	<< endl;
gsell's avatar
gsell committed
171 172 173 174 175 176 177 178 179 180 181
}

void ParallelSliceTracker::applyEntranceFringe(double angle, double curve,
        const BMultipoleField &field, double scale) {
}


void ParallelSliceTracker::applyExitFringe(double angle, double curve,
        const BMultipoleField &field, double scale) {
}

182

gsell's avatar
gsell committed
183 184
void ParallelSliceTracker::execute() {

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

189
    currentSimulationTime_m = itsBunch_m->getT();
gsell's avatar
gsell committed
190

191 192 193 194
    setTime();
    setLastStep();

    prepareSections();
195
    handleAutoPhasing();
gsell's avatar
gsell committed
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 232 233 234 235 236 237 238 239 240 241
    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
242
    Vector_t rmin, rmax;
243 244 245 246 247 248 249 250 251 252 253 254 255 256 257
    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
258

259
    setTime();
gsell's avatar
gsell committed
260

261 262 263 264 265 266 267
    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
268

269 270 271 272 273
    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
274

275
    // prepareEmission();
276

277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293
    *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) {

294
            globalEOL_m = false;
295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312

            switchElements();
            computeExternalFields(oth);

            //reduce(&globalEOL_m, &globalEOL_m, OpBitwiseOrAssign());
            //reduce(&globalEOL_m, &globalEOL_m + 1, &globalEOL_m, OpBitwiseAndAssign());
            MPI_Allreduce(MPI_IN_PLACE, &globalEOL_m, 1, MPI_INT, MPI_LAND, Ippl::getComm());

            computeSpaceChargeFields();
            timeIntegration();

            if (step % 1000 + 1 == 1000) {
                msg << level1;
            } else if (step % 100 + 1 == 100) {
                msg << level2;
            } else {
                msg << level3;
            }
313 314 315 316 317

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

318 319
            currentSimulationTime_m += itsBunch_m->getdT();
            itsBunch_m->setT(currentSimulationTime_m);
320

321 322
            dumpStats(step);

323 324 325 326 327 328 329 330 331
            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;
332
        }
333

334
        if (globalEOL_m)
335 336 337 338 339
            break;

        dtAllTracks_m.pop();
        localTrackSteps_m.pop();
        zStop_m.pop();
gsell's avatar
gsell committed
340 341
    }

342 343
    OpalData::getInstance()->setLastStep(step);
    writeLastStepPhaseSpace(step, itsBunch_m->get_sPos());
344
    itsOpalBeamline_m.switchElementsOff();
345 346 347 348 349 350 351 352 353 354 355 356
    msg << "done executing ParallelSliceTracker" << endl;
}


void ParallelSliceTracker::timeIntegration() {

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


357
void ParallelSliceTracker::handleAutoPhasing() {
358

359
    if (Options::autoPhase == 0) return;
gsell's avatar
gsell committed
360

361
    if(!OpalData::getInstance()->inRestartRun()) {
362
        itsDataSink_m->storeCavityInformation();
gsell's avatar
gsell committed
363 364
    }

365 366 367 368 369 370
    auto it = OpalData::getInstance()->getFirstMaxPhases();
    auto end = OpalData::getInstance()->getLastMaxPhases();
    for (; it != end; ++ it) {
        updateRFElement((*it).first, (*it).second);
    }
    printRFPhases();
371
}
gsell's avatar
gsell committed
372 373


374
void ParallelSliceTracker::prepareSections() {
gsell's avatar
gsell committed
375

376
    itsBeamline_m.accept(*this);
377
    itsOpalBeamline_m.prepareSections();
gsell's avatar
gsell committed
378

379 380 381
    itsOpalBeamline_m.compute3DLattice();
    itsOpalBeamline_m.save3DLattice();
    itsOpalBeamline_m.save3DInput();
382
}
gsell's avatar
gsell committed
383 384


385
void ParallelSliceTracker::computeExternalFields(OrbitThreader &oth) {
gsell's avatar
gsell committed
386

387
    IpplTimings::startTimer(timeFieldEvaluation_m);
gsell's avatar
gsell committed
388

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

391 392 393 394 395 396 397 398 399 400 401 402
    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
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 431 432 433 434 435 436 437 438 439 440
    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
441

442 443
            (*it)->addKR(i, time, KRThis);
            (*it)->addKT(i, time, KTThis);
gsell's avatar
gsell committed
444

445 446
            externalE += fromLocal.rotateTo(externalEThis);
            externalB += fromLocal.rotateTo(externalBThis);
gsell's avatar
gsell committed
447

448 449
            KR += fromLocal.rotateTo(KRThis);
            KT += fromLocal.rotateTo(KTThis);
gsell's avatar
gsell committed
450

451 452 453
            itsBunch_m->setR(i, pos);
            itsBunch_m->setP(i, mom);
        }
gsell's avatar
gsell committed
454

455 456
        itsBunch_m->setExternalFields(i, externalE, externalB, KR, KT);
    }
gsell's avatar
gsell committed
457

458 459 460 461 462 463 464 465
    IpplTimings::stopTimer(timeFieldEvaluation_m);
}


void ParallelSliceTracker::computeSpaceChargeFields() {

    itsBunch_m->computeSpaceCharge();
}
gsell's avatar
gsell committed
466 467


468
void ParallelSliceTracker::dumpStats(long long step) {
gsell's avatar
gsell committed
469

snuverink_j's avatar
snuverink_j committed
470 471 472
    double sposRef = itsBunch_m->get_sPos();
    if (step != 0 && (step % Options::psDumpFreq == 0 || step % Options::statDumpFreq == 0))
        writePhaseSpace(step, sposRef);
473
}
gsell's avatar
gsell committed
474

475 476 477

void ParallelSliceTracker::switchElements(double scaleMargin) {

snuverink_j's avatar
snuverink_j committed
478
    double margin = 1.0;
479

snuverink_j's avatar
snuverink_j committed
480
    currentSimulationTime_m = itsBunch_m->getT();
481
    itsOpalBeamline_m.switchElements(itsBunch_m->zTail() - margin,
snuverink_j's avatar
snuverink_j committed
482 483
                                      itsBunch_m->zHead() + margin,
                                      itsBunch_m->Eavg() * 1e-6);
484 485 486 487 488
}


void ParallelSliceTracker::setLastStep() {

489
    // unsigned long long step = 0;
490

491
    // if (OpalData::getInstance()->inRestartRun()) {
492

493 494 495 496
    //     int prevDumpFreq = OpalData::getInstance()->getRestartDumpFreq();
    //     step = OpalData::getInstance()->getRestartStep() * prevDumpFreq + 1;
    //     maxSteps_m += step;
    // } else {
497

498 499 500
    //     step = OpalData::getInstance()->getLastStep() + 1;
    //     maxSteps_m += step;
    // }
gsell's avatar
gsell committed
501

502
    // OpalData::getInstance()->setLastStep(step);
503
}
gsell's avatar
gsell committed
504 505


506
void ParallelSliceTracker::setTime() {
gsell's avatar
gsell committed
507

508 509 510 511 512 513 514 515
    if (OpalData::getInstance()->inRestartRun())
        currentSimulationTime_m = itsBunch_m->getT();
    else
        currentSimulationTime_m = itsBunch_m->getT();
}


bool ParallelSliceTracker::hasEndOfLineReached() {
516 517
    reduce(&globalEOL_m, &globalEOL_m + 1, &globalEOL_m, OpBitwiseAndAssign());
    return globalEOL_m;
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 667 668 669 670 671 672 673 674 675 676

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();
    }
}