ParallelCyclotronTracker.cpp 125 KB
Newer Older
gsell's avatar
gsell committed
1
//
2 3
// Class ParallelCyclotronTracker
//   Tracker for OPAL-Cycl
gsell's avatar
gsell committed
4
//
5 6 7 8 9
// Copyright (c) 2007 - 2014, Jianjun Yang, Andreas Adelmann and Matthias Toggweiler,
//                            Paul Scherrer Institut, Villigen PSI, Switzerland
// Copyright (c) 2014,        Daniel Winklehner, MIT, Cambridge, MA, USA
// Copyright (c) 2012 - 2020, Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved
gsell's avatar
gsell committed
10
//
11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
// Implemented as part of the PhD thesis
// "Beam dynamics in high intensity cyclotrons including neighboring bunch effects"
// and the paper
// "Beam dynamics in high intensity cyclotrons including neighboring bunch effects:
// Model, implementation, and application"
// (https://journals.aps.org/prab/pdf/10.1103/PhysRevSTAB.13.064201)
//
// 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/>.
gsell's avatar
gsell committed
27
//
kraus's avatar
kraus committed
28
#include "Algorithms/ParallelCyclotronTracker.h"
29

gsell's avatar
gsell committed
30
#include <fstream>
31 32
#include <iostream>
#include <limits>
gsell's avatar
gsell committed
33
#include <vector>
frey_m's avatar
frey_m committed
34
#include <numeric>
35
#include <cmath>
36 37

#include "AbstractObjects/Element.h"
38
#include "AbstractObjects/OpalData.h"
gsell's avatar
gsell committed
39

40
#include "AbsBeamline/BeamStripping.h"
41
#include "AbsBeamline/CCollimator.h"
gsell's avatar
gsell committed
42 43
#include "AbsBeamline/Corrector.h"
#include "AbsBeamline/Cyclotron.h"
adelmann's avatar
adelmann committed
44
#include "AbsBeamline/Degrader.h"
gsell's avatar
gsell committed
45 46
#include "AbsBeamline/Drift.h"
#include "AbsBeamline/ElementBase.h"
47
#include "AbsBeamline/Offset.h"
gsell's avatar
gsell committed
48 49 50
#include "AbsBeamline/Marker.h"
#include "AbsBeamline/Monitor.h"
#include "AbsBeamline/Multipole.h"
ext-rogers_c's avatar
ext-rogers_c committed
51
#include "AbsBeamline/MultipoleT.h"
52 53 54 55
#include "AbsBeamline/MultipoleTBase.h"
#include "AbsBeamline/MultipoleTStraight.h"
#include "AbsBeamline/MultipoleTCurvedConstRadius.h"
#include "AbsBeamline/MultipoleTCurvedVarRadius.h"
56
#include "AbsBeamline/PluginElement.h"
gsell's avatar
gsell committed
57 58
#include "AbsBeamline/Probe.h"
#include "AbsBeamline/RBend.h"
59
#include "AbsBeamline/Ring.h"
gsell's avatar
gsell committed
60 61
#include "AbsBeamline/RFCavity.h"
#include "AbsBeamline/SBend.h"
62
#include "AbsBeamline/SBend3D.h"
ext-rogers_c's avatar
ext-rogers_c committed
63
#include "AbsBeamline/ScalingFFAMagnet.h"
gsell's avatar
gsell committed
64 65 66
#include "AbsBeamline/Septum.h"
#include "AbsBeamline/Solenoid.h"
#include "AbsBeamline/Stripper.h"
67
#include "AbsBeamline/VariableRFCavity.h"
68
#include "AbsBeamline/VariableRFCavityFringeField.h"
69
#include "AbsBeamline/VerticalFFAMagnet.h"
70

71 72
#include "Algorithms/Ctunes.h"
#include "Algorithms/PolynomialTimeDependence.h"
gsell's avatar
gsell committed
73 74

#include "Beamlines/Beamline.h"
75
#include "Beamlines/FlaggedBeamline.h"
gsell's avatar
gsell committed
76

77
#include "Elements/OpalBeamline.h"
gsell's avatar
gsell committed
78 79 80 81

#include "Physics/Physics.h"

#include "Utilities/OpalException.h"
82
#include "Utilities/Options.h"
gsell's avatar
gsell committed
83

84
#include "BasicActions/DumpFields.h"
85
#include "BasicActions/DumpEMFields.h"
86

87
#include "Structure/BoundaryGeometry.h"
88
#include "Structure/DataSink.h"
Pedro Calvo Portela's avatar
Pedro Calvo Portela committed
89
#include "Structure/LossDataSink.h"
frey_m's avatar
frey_m committed
90

91

92
constexpr double c_mmtns = Physics::c * 1.0e-6; // m/s --> mm/ns
93

94 95 96
Vector_t const ParallelCyclotronTracker::xaxis = Vector_t(1.0, 0.0, 0.0);
Vector_t const ParallelCyclotronTracker::yaxis = Vector_t(0.0, 1.0, 0.0);
Vector_t const ParallelCyclotronTracker::zaxis = Vector_t(0.0, 0.0, 1.0);
gsell's avatar
gsell committed
97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112

extern Inform *gmsg;

/**
 * Constructor ParallelCyclotronTracker
 *
 * @param beamline
 * @param bunch
 * @param ds
 * @param reference
 * @param revBeam
 * @param revTrack
 * @param maxSTEPS
 * @param timeIntegrator
 */
ParallelCyclotronTracker::ParallelCyclotronTracker(const Beamline &beamline,
frey_m's avatar
frey_m committed
113
                                                   PartBunchBase<double, 3> *bunch,
gsell's avatar
gsell committed
114 115 116
                                                   DataSink &ds,
                                                   const PartData &reference,
                                                   bool revBeam, bool revTrack,
frey_m's avatar
frey_m committed
117
                                                   int maxSTEPS, int timeIntegrator,
118 119 120 121 122
                                                   const int& numBunch,
                                                   const double& mbEta,
                                                   const double& mbPara,
                                                   const std::string& mbMode,
                                                   const std::string& mbBinning)
frey_m's avatar
frey_m committed
123 124 125 126 127 128 129 130 131 132
    : Tracker(beamline, bunch, reference, revBeam, revTrack)
    , bgf_m(nullptr)
    , maxSteps_m(maxSTEPS)
    , lastDumpedStep_m(0)
    , myNode_m(Ippl::myNode())
    , initialLocalNum_m(bunch->getLocalNum())
    , initialTotalNum_m(bunch->getTotalNum())
    , opalRing_m(nullptr)
    , itsStepper_mp(nullptr)
{
gsell's avatar
gsell committed
133 134 135
    itsBeamline = dynamic_cast<Beamline *>(beamline.clone());
    itsDataSink = &ds;

136 137 138
    if ( numBunch > 1 ) {
        mbHandler_m = std::unique_ptr<MultiBunchHandler>(
            new MultiBunchHandler(bunch, numBunch, mbEta,
frey_m's avatar
frey_m committed
139
                                  mbPara, mbMode, mbBinning)
140 141 142
        );
    }

gsell's avatar
gsell committed
143 144 145 146
    IntegrationTimer_m = IpplTimings::getTimer("Integration");
    TransformTimer_m   = IpplTimings::getTimer("Frametransform");
    DumpTimer_m        = IpplTimings::getTimer("Dump");
    BinRepartTimer_m   = IpplTimings::getTimer("Binaryrepart");
147 148
    PluginElemTimer_m  = IpplTimings::getTimer("PluginElements");
    DelParticleTimer_m = IpplTimings::getTimer("DeleteParticles");
149

150 151 152 153 154 155 156 157 158
    // FIXME Change track command
    if ( initialTotalNum_m == 1 ) {
        mode_m = MODE::SINGLE;
    } else if ( initialTotalNum_m == 2 ) {
        mode_m = MODE::SEO;
    } else if ( initialTotalNum_m > 2 ) {
        mode_m = MODE::BUNCH;
    } else
        mode_m = MODE::UNDEFINED;
159

160 161 162 163 164 165 166 167
    if ( timeIntegrator == 0 ) {
        stepper_m = stepper::INTEGRATOR::RK4;
    } else if ( timeIntegrator == 1) {
        stepper_m = stepper::INTEGRATOR::LF2;
    } else if ( timeIntegrator == 2) {
        stepper_m = stepper::INTEGRATOR::MTS;
    } else
        stepper_m = stepper::INTEGRATOR::UNDEFINED;
gsell's avatar
gsell committed
168 169 170 171 172 173 174
}

/**
 * Destructor ParallelCyclotronTracker
 *
 */
ParallelCyclotronTracker::~ParallelCyclotronTracker() {
Pedro Calvo Portela's avatar
Pedro Calvo Portela committed
175 176
    if(bgf_m)
        lossDs_m->save();
177 178
    for(Component* component : myElements) {
        delete(component);
gsell's avatar
gsell committed
179
    }
180 181
    for(auto fd : FieldDimensions) {
        delete(fd);
gsell's avatar
gsell committed
182 183
    }
    delete itsBeamline;
184
    // delete opalRing_m;
gsell's avatar
gsell committed
185 186
}

frey_m's avatar
frey_m committed
187

188
void ParallelCyclotronTracker::bgf_main_collision_test() {
189
    if(!bgf_m) return;
190

191
    Inform msg("bgf_main_collision_test ");
192

193
    /**
194
     *Here we check if a particle is outside the domain, flag it for deletion
195
     */
196

197
    Vector_t intecoords = 0.0;
198

199
    // This has to match the dT in the rk4 pusher
frey_m's avatar
frey_m committed
200
    double dtime = itsBunch_m->getdT() * getHarmonicNumber();
201

202
    int triId = 0;
frey_m's avatar
frey_m committed
203
    for(size_t i = 0; i < itsBunch_m->getLocalNum(); i++) {
204 205
        int res = bgf_m->partInside(itsBunch_m->R[i], itsBunch_m->P[i],
                                    dtime, intecoords, triId);
206
        if(res >= 0) {
ext-calvo_p's avatar
ext-calvo_p committed
207
            lossDs_m->addParticle(itsBunch_m->R[i], itsBunch_m->P[i],
Pedro Calvo Portela's avatar
Pedro Calvo Portela committed
208 209
                                  itsBunch_m->ID[i], itsBunch_m->getT()*1e9,
                                  turnnumber_m, itsBunch_m->bunchNum[i]);
frey_m's avatar
frey_m committed
210
            itsBunch_m->Bin[i] = -1;
211
            Inform gmsgALL("OPAL", INFORM_ALL_NODES);
212 213
            gmsgALL << level4 << "* Particle " << itsBunch_m->ID[i]
                    << " lost on boundary geometry" << endl;
214
        }
215
    }
216 217
}

218
// only used for dumping into stat file
frey_m's avatar
frey_m committed
219 220
void ParallelCyclotronTracker::dumpAngle(const double& theta,
                                         double& prevAzimuth,
frey_m's avatar
frey_m committed
221 222
                                         double& azimuth,
                                         const short& bunchNr)
frey_m's avatar
frey_m committed
223 224
{
    if ( prevAzimuth < 0.0 ) { // only at first occurrence
frey_m's avatar
frey_m committed
225 226
        double plus = 0.0;
        if ( OpalData::getInstance()->inRestartRun() ) {
frey_m's avatar
frey_m committed
227
            plus = 360.0 * (turnnumber_m - bunchNr);
frey_m's avatar
frey_m committed
228
        }
frey_m's avatar
frey_m committed
229
        azimuth = theta + plus;
230
    } else {
frey_m's avatar
frey_m committed
231
        double dtheta = theta - prevAzimuth;
232 233 234 235 236 237
        if ( dtheta < 0 ) {
            dtheta += 360.0;
        }
        if ( dtheta > 180 ) { // rotating clockwise, reduce angle
            dtheta -= 360;
        }
frey_m's avatar
frey_m committed
238
        azimuth += dtheta;
239
    }
frey_m's avatar
frey_m committed
240
    prevAzimuth = theta;
241 242 243
}


frey_m's avatar
frey_m committed
244 245 246 247 248 249
double ParallelCyclotronTracker::computeRadius(const Vector_t& meanR) const {
    // New OPAL 2.0: m --> mm
    return 1000.0 * std::sqrt(meanR(0) * meanR(0) + meanR(1) * meanR(1));
}


frey_m's avatar
frey_m committed
250 251
void ParallelCyclotronTracker::computePathLengthUpdate(std::vector<double>& dl,
                                                       const double& dt)
252
{
frey_m's avatar
frey_m committed
253 254
    // the last element in dotP is the dot-product over all particles
    std::vector<double> dotP(dl.size());
255
    if ( Options::psDumpFrame == Options::BUNCH_MEAN || isMultiBunch()) {
256

frey_m's avatar
frey_m committed
257 258
        for(unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) {
            dotP[itsBunch_m->bunchNum[i]] += dot(itsBunch_m->P[i], itsBunch_m->P[i]);
259 260
        }

frey_m's avatar
frey_m committed
261
        allreduce(dotP.data(), dotP.size(), std::plus<double>());
262

frey_m's avatar
frey_m committed
263 264 265
        // dot-product over all particles
        double sum = std::accumulate(dotP.begin(), dotP.end() - 1, 0);
        dotP.back() = sum / double(itsBunch_m->getTotalNum());
frey_m's avatar
frey_m committed
266

frey_m's avatar
frey_m committed
267 268 269 270
        // bunch specific --> multi-bunches only
        for (short b = 0; b < (short)dotP.size() - 1; ++b) {
            dotP[b] /= double(itsBunch_m->getTotalNumPerBunch(b));
        }
271 272

    } else if ( itsBunch_m->getLocalNum() == 0 ) {
273
        // here we are in Options::GLOBAL mode
frey_m's avatar
frey_m committed
274
        dotP[0] = 0.0;
275
    } else {
276
        // here we are in Options::GLOBAL mode
frey_m's avatar
frey_m committed
277
        dotP[0] = dot(itsBunch_m->P[0], itsBunch_m->P[0]);
278 279
    }

frey_m's avatar
frey_m committed
280 281 282 283 284
    for (size_t i = 0; i < dotP.size(); ++i) {
        double const gamma = std::sqrt(1.0 + dotP[i]);
        double const beta  = std::sqrt(dotP[i]) / gamma;
        dl[i] = c_mmtns * dt * 1.0e-3 * beta; // unit: m
    }
285 286 287
}


288

gsell's avatar
gsell committed
289 290 291 292 293
/**
 *
 *
 * @param fn Base file name
 */
294
void ParallelCyclotronTracker::openFiles(size_t numFiles, std::string SfileName) {
gsell's avatar
gsell committed
295

296
    for (unsigned int i=0; i<numFiles; i++) {
297 298 299 300 301 302 303 304
        std::string SfileName2 = SfileName;
        if (i<=2) {
            SfileName2 += std::string("-Angle" + std::to_string(int(i)) + ".dat");
        }
        else {
            // for single Particle Mode, output after each turn, to define matched initial phase ellipse.
            SfileName2 += std::string("-afterEachTurn.dat");
        }
gsell's avatar
gsell committed
305

306 307 308 309 310 311
        outfTheta_m.emplace_back(new std::ofstream(SfileName2.c_str()));
        outfTheta_m.back()->precision(8);
        outfTheta_m.back()->setf(std::ios::scientific, std::ios::floatfield);
        *outfTheta_m.back() << "# r [mm]        beta_r*gamma       "
                            << "theta [deg]     beta_theta*gamma        "
                            << "z [mm]          beta_z*gamma" << std::endl;
312
    }
gsell's avatar
gsell committed
313 314 315 316 317 318 319 320
}

/**
 * Close all files related to
 * special output in the Cyclotron
 * mode.
 */
void ParallelCyclotronTracker::closeFiles() {
321
    for (auto & file : outfTheta_m) {
322
        file->close();
323
    }
gsell's avatar
gsell committed
324 325
}

326
/**
327 328 329
 *
 * @param ring
 */
kraus's avatar
kraus committed
330
void ParallelCyclotronTracker::visitRing(const Ring &ring) {
Daniel Winklehner's avatar
Daniel Winklehner committed
331

332
    *gmsg << "* ----------------------------- Adding Ring ------------------------------ *" << endl;
Daniel Winklehner's avatar
Daniel Winklehner committed
333

snuverink_j's avatar
snuverink_j committed
334
    delete opalRing_m;
Daniel Winklehner's avatar
Daniel Winklehner committed
335

kraus's avatar
kraus committed
336
    opalRing_m = dynamic_cast<Ring*>(ring.clone());
Daniel Winklehner's avatar
Daniel Winklehner committed
337

338
    myElements.push_back(opalRing_m);
Daniel Winklehner's avatar
Daniel Winklehner committed
339

frey_m's avatar
frey_m committed
340
    opalRing_m->initialise(itsBunch_m);
341 342 343 344

    referenceR = opalRing_m->getBeamRInit();
    referencePr = opalRing_m->getBeamPRInit();
    referenceTheta = opalRing_m->getBeamPhiInit();
Daniel Winklehner's avatar
Daniel Winklehner committed
345

346
    if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
kraus's avatar
kraus committed
347
        throw OpalException("Error in ParallelCyclotronTracker::visitRing",
348 349
                            "PHIINIT is out of [-180, 180)!");
    }
Daniel Winklehner's avatar
Daniel Winklehner committed
350 351

    referenceZ = 0.0;
352
    referencePz = 0.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
353 354

    referencePtot = itsReference.getGamma() * itsReference.getBeta();
355
    referencePt = std::sqrt(referencePtot * referencePtot - referencePr * referencePr);
Daniel Winklehner's avatar
Daniel Winklehner committed
356

357 358
    if(referencePtot < 0.0)
        referencePt *= -1.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
359

360 361
    sinRefTheta_m = std::sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = std::cos(referenceTheta * Physics::deg2rad);
gsell's avatar
gsell committed
362

snuverink_j's avatar
snuverink_j committed
363
    double BcParameter[8] = {}; // zero initialise array
Daniel Winklehner's avatar
Daniel Winklehner committed
364

365
    buildupFieldList(BcParameter, ElementBase::RING, opalRing_m);
366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381

    // Finally print some diagnostic
    *gmsg << "* Initial beam radius = " << referenceR << " [mm] " << endl;
    *gmsg << "* Initial gamma = " << itsReference.getGamma() << endl;
    *gmsg << "* Initial beta = " << itsReference.getBeta() << endl;
    *gmsg << "* Total reference momentum   = " << referencePtot * 1000.0
          << " [MCU]" << endl;
    *gmsg << "* Reference azimuthal momentum  = " << referencePt * 1000.0
          << " [MCU]" << endl;
    *gmsg << "* Reference radial momentum     = " << referencePr * 1000.0
          << " [MCU]" << endl;
    *gmsg << "* " << opalRing_m->getSymmetry() << " fold field symmetry "
          << endl;
    *gmsg << "* Harmonic number h= " << opalRing_m->getHarmonicNumber() << " "
          << endl;
}
gsell's avatar
gsell committed
382 383 384 385 386 387 388 389

/**
 *
 *
 * @param cycl
 */
void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) {

390
    *gmsg << "* -------------------------- Adding Cyclotron ---------------------------- *" << endl;
gsell's avatar
gsell committed
391

392 393
    cycl_m = dynamic_cast<Cyclotron *>(cycl.clone());
    myElements.push_back(cycl_m);
394

395
    // Is this a Spiral Inflector Simulation? If yes, we'll give the user some
396
    // useful information
397
    spiral_flag = cycl_m->getSpiralFlag();
398 399 400 401

    if(spiral_flag) {

        *gmsg << endl << "* This is a Spiral Inflector Simulation! This means the following:" << endl;
402
        *gmsg         << "* 1.) It is up to the user to provide appropriate geometry, electric and magnetic fields!" << endl;
403 404
        *gmsg         << "*     (Use BANDRF type cyclotron and use RFMAPFN to load both magnetic" << endl;
        *gmsg         << "*     and electric fields, setting SUPERPOSE to an array of TRUE values.)" << endl;
snuverink_j's avatar
snuverink_j committed
405
        *gmsg         << "* 2.) For high currents it is strongly recommended to use the SAAMG fieldsolver," << endl;
406 407
        *gmsg         << "*     FFT does not give the correct results (boundary conditions are missing)." << endl;
        *gmsg         << "* 3.) The whole geometry will be meshed and used for the fieldsolver." << endl;
408 409
        *gmsg         << "*     There will be no transformations of the bunch into a local frame und consequently," << endl;
        *gmsg         << "*     the problem will be treated non-relativistically!" << endl;
410 411 412
        *gmsg         << "*     (This is not an issue for spiral inflectors as they are typically < 100 keV/amu.)" << endl;
        *gmsg << endl << "* Note: For now, multi-bunch mode (MBM) needs to be de-activated for spiral inflector" << endl;
        *gmsg         << "* and space charge needs to be solved every time-step. numBunch_m and scSolveFreq are reset." << endl;
frey_m's avatar
frey_m committed
413 414 415
        if (isMultiBunch()) {
            mbHandler_m = nullptr;
        }
416 417
    }

418
    // Fresh run (no restart):
419
    if(!OpalData::getInstance()->inRestartRun()) {
420

421
        // Get reference values from cyclotron element
422
        // For now, these are still stored in mm. should be the only ones. -DW
423 424 425 426 427
        referenceR     = cycl_m->getRinit();
        referenceTheta = cycl_m->getPHIinit();
        referenceZ     = cycl_m->getZinit();
        referencePr    = cycl_m->getPRinit();
        referencePz    = cycl_m->getPZinit();
428

429
        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
430 431
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                "PHIINIT is out of [-180, 180)!");
432 433 434 435 436
        }

        referencePtot =  itsReference.getGamma() * itsReference.getBeta();

        // Calculate reference azimuthal (tangential) momentum from total-, z- and radial momentum:
437
        float insqrt = referencePtot * referencePtot - \
438
            referencePr * referencePr - referencePz * referencePz;
439 440 441 442 443

        if(insqrt < 0) {

            if(insqrt > -1.0e-10) {

444
                referencePt = 0.0;
445 446 447

            } else {

448
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
449
                                    "Pt imaginary!");
450 451 452 453
            }

        } else {

454
            referencePt = std::sqrt(insqrt);
455 456
        }

457
        if(referencePtot < 0.0)
Daniel Winklehner's avatar
Daniel Winklehner committed
458
            referencePt *= -1.0;
459 460
        // End calculate referencePt

461
        // Restart a run:
462 463
    } else {

464
        // If the user wants to save the restarted run in local frame,
465
        // make sure the previous h5 file was local too
466
      if (Options::psDumpFrame != Options::GLOBAL) {
winklehner_d's avatar
-DW  
winklehner_d committed
467
        if (!previousH5Local) {
468 469
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                    "You are trying a local restart from a global h5 file!");
winklehner_d's avatar
-DW  
winklehner_d committed
470
        }
471 472
            // Else, if the user wants to save the restarted run in global frame,
            // make sure the previous h5 file was global too
winklehner_d's avatar
-DW  
winklehner_d committed
473 474
      } else {
        if (previousH5Local) {
475 476
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                    "You are trying a global restart from a local h5 file!");
477
            }
478
        }
479

480
        // Adjust some of the reference variables from the h5 file
481 482
        referencePhi *= Physics::deg2rad;
        referencePsi *= Physics::deg2rad;
483
        referencePtot = bega;
484
        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
485
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
486 487
                                "PHIINIT is out of [-180, 180)!");
        }
488 489
    }

490 491
    sinRefTheta_m = std::sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = std::cos(referenceTheta * Physics::deg2rad);
492

493
    *gmsg << endl;
adelmann's avatar
adelmann committed
494
    *gmsg << "* Bunch global starting position:" << endl;
495 496
    *gmsg << "* RINIT = " << referenceR  << " [mm]" << endl;
    *gmsg << "* PHIINIT = " << referenceTheta << " [deg]" << endl;
497
    *gmsg << "* ZINIT = " << referenceZ << " [mm]" << endl;
498
    *gmsg << endl;
adelmann's avatar
adelmann committed
499
    *gmsg << "* Bunch global starting momenta:" << endl;
gsell's avatar
gsell committed
500 501
    *gmsg << "* Initial gamma = " << itsReference.getGamma() << endl;
    *gmsg << "* Initial beta = " << itsReference.getBeta() << endl;
502
    *gmsg << "* Reference total momentum (beta * gamma) = " << referencePtot * 1000.0 << " [MCU]" << endl;
503 504 505
    *gmsg << "* Reference azimuthal momentum (Pt) = " << referencePt * 1000.0 << " [MCU]" << endl;
    *gmsg << "* Reference radial momentum (Pr) = " << referencePr * 1000.0 << " [MCU]" << endl;
    *gmsg << "* Reference axial momentum (Pz) = " << referencePz * 1000.0 << " [MCU]" << endl;
506
    *gmsg << endl;
adelmann's avatar
adelmann committed
507

508
    double sym = cycl_m->getSymmetry();
snuverink_j's avatar
snuverink_j committed
509
    *gmsg << "* " << sym << "-fold field symmetry " << endl;
gsell's avatar
gsell committed
510

511
    // ckr: this just returned the default value as defined in Component.h
512
    // double rff = cycl_m->getRfFrequ(0);
513
    // *gmsg << "* Rf frequency= " << rff << " [MHz]" << endl;
gsell's avatar
gsell committed
514

515
    std::string fmfn = cycl_m->getFieldMapFN();
516
    *gmsg << "* Field map file name = " << fmfn << " " << endl;
gsell's avatar
gsell committed
517

518
    std::string type = cycl_m->getCyclotronType();
519
    *gmsg << "* Type of cyclotron = " << type << " " << endl;
520

521 522
    double rmin = cycl_m->getMinR();
    double rmax = cycl_m->getMaxR();
adelmann's avatar
adelmann committed
523
    *gmsg << "* Radial aperture = " << rmin << " ... " << rmax<<" [m] "<< endl;
524

525 526
    double zmin = cycl_m->getMinZ();
    double zmax = cycl_m->getMaxZ();
adelmann's avatar
adelmann committed
527
    *gmsg << "* Vertical aperture = " << zmin << " ... " << zmax<<" [m]"<< endl;
gsell's avatar
gsell committed
528

529 530
    double h = cycl_m->getCyclHarm();
    *gmsg << "* Number of trimcoils = " << cycl_m->getNumberOfTrimcoils() << endl;
531
    *gmsg << "* Harmonic number h = " << h << " " << endl;
532

533 534 535 536 537 538 539 540 541
    if (type == std::string("BANDRF")) {
        double escale = cycl_m->getEScale(0);
        *gmsg << "* RF field scale factor = " << escale << endl;
        double rfphi= cycl_m->getRfPhi(0);
        *gmsg << "* RF inital phase = " << rfphi * Physics::rad2deg << " [deg]" << endl;
        bool superpose = cycl_m->getSuperpose(0);
        *gmsg << std::boolalpha << "* Superpose electric field maps -> " << superpose << endl;
    }

gsell's avatar
gsell committed
542 543 544 545 546 547 548
    /**
     * To ease the initialise() function, set a integral parameter fieldflag internally.
     * Its value is  by the option "TYPE" of the element  "CYCLOTRON"
     * fieldflag = 1, readin PSI format measured field file (default)
     * fieldflag = 2, readin carbon cyclotron field file created by Jianjun Yang, TYPE=CARBONCYCL
     * fieldflag = 3, readin ANSYS format file for CYCIAE-100 created by Jianjun Yang, TYPE=CYCIAE
     * fieldflag = 4, readin AVFEQ format file for Riken cyclotrons
ext-rogers_c's avatar
ext-rogers_c committed
549
     * fieldflag = 5, readin FFA format file for MSU/FNAL FFA
gsell's avatar
gsell committed
550
     * fieldflag = 6, readin both median plane B field map and 3D E field map of RF cavity for compact cyclotron
551
     * fieldflag = 7, read in fields for Daniel's synchrocyclotron simulations
gsell's avatar
gsell committed
552
     */
553
    int fieldflag = cycl_m->getFieldFlag(type);
gsell's avatar
gsell committed
554

555
    // Read in cyclotron field maps (midplane + 3D fields if desired).
556
    cycl_m->initialise(itsBunch_m, fieldflag, cycl_m->getBScale());
gsell's avatar
gsell committed
557

snuverink_j's avatar
snuverink_j committed
558
    double BcParameter[8] = {};
559

560 561
    BcParameter[0] = 0.001 * cycl_m->getRmin();
    BcParameter[1] = 0.001 * cycl_m->getRmax();
gsell's avatar
gsell committed
562

563
    // Store inner radius and outer radius of cyclotron field map in the list
564
    buildupFieldList(BcParameter, ElementBase::CYCLOTRON, cycl_m);
gsell's avatar
gsell committed
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

void ParallelCyclotronTracker::visitBeamStripping(const BeamStripping &bstp) {
    *gmsg << "* ------------------------------ Beam Stripping ------------------------------" << endl;

    BeamStripping* elptr = dynamic_cast<BeamStripping *>(bstp.clone());
    myElements.push_back(elptr);

    double BcParameter[8] = {};

    if(elptr->getPressureMapFN() == "") {
        double pressure = elptr->getPressure();
        *gmsg << "* Pressure     = " << pressure << " [mbar]" << endl;
        BcParameter[0] = pressure;
    }

    double temperature = elptr->getTemperature();
    *gmsg << "* Temperature  = " << temperature << " [K]" << endl;

    std::string gas = elptr->getResidualGas();
    *gmsg << "* Residual gas = " << gas << endl;

    bool stop = elptr->getStop();
    *gmsg << std::boolalpha << "* Particles stripped will be deleted after interaction -> " << stop << endl;

    elptr->initialise(itsBunch_m, elptr->getPScale());
    
    BcParameter[1] = temperature;
    BcParameter[2] = stop;

    buildupFieldList(BcParameter, ElementBase::BEAMSTRIPPING, elptr);

}

gsell's avatar
gsell committed
600 601 602 603 604
/**
 *
 *
 * @param coll
 */
605
void ParallelCyclotronTracker::visitCCollimator(const CCollimator &coll) {
gsell's avatar
gsell committed
606

607
    *gmsg << "* ------------------------------ Collimator ------------------------------" << endl;
gsell's avatar
gsell committed
608

609
    CCollimator* elptr = dynamic_cast<CCollimator *>(coll.clone());
610
    myElements.push_back(elptr);
gsell's avatar
gsell committed
611

612
    double xstart = elptr->getXStart();
613
    *gmsg << "* Xstart  = " << xstart << " [m]" << endl;
gsell's avatar
gsell committed
614

615
    double xend = elptr->getXEnd();
616
    *gmsg << "* Xend    = " << xend << " [m]" << endl;
gsell's avatar
gsell committed
617

618
    double ystart = elptr->getYStart();
619
    *gmsg << "* Ystart  = " << ystart << " [m]" << endl;
gsell's avatar
gsell committed
620

621
    double yend = elptr->getYEnd();
622
    *gmsg << "* Yend    = " << yend << " [m]" << endl;
gsell's avatar
gsell committed
623

624
    double zstart = elptr->getZStart();
625
    *gmsg << "* Zstart  = " << zstart << " [m]" << endl;
626 627

    double zend = elptr->getZEnd();
628
    *gmsg << "* Zend    = " << zend << " [m]" << endl;
629

630
    double width = elptr->getWidth();
631
    *gmsg << "* Width   = " << width << " [m]" << endl;
gsell's avatar
gsell committed
632

frey_m's avatar
frey_m committed
633
    elptr->initialise(itsBunch_m);
gsell's avatar
gsell committed
634

snuverink_j's avatar
snuverink_j committed
635
    double BcParameter[8] = {};
636

637 638 639 640 641
    BcParameter[0] = xstart;
    BcParameter[1] = xend;
    BcParameter[2] = ystart;
    BcParameter[3] = yend;
    BcParameter[4] = width;
642

643
    buildupFieldList(BcParameter, ElementBase::CCOLLIMATOR, elptr);
gsell's avatar
gsell committed
644 645 646 647 648 649 650 651 652 653 654 655
}

/**
 *
 *
 * @param corr
 */
void ParallelCyclotronTracker::visitCorrector(const Corrector &corr) {
    *gmsg << "In Corrector; L= " << corr.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Corrector *>(corr.clone()));
}

adelmann's avatar
adelmann committed
656 657 658 659 660 661 662 663 664 665 666 667
/**
 *
 *
 * @param degrader
 */
void ParallelCyclotronTracker::visitDegrader(const Degrader &deg) {
    *gmsg << "In Degrader; L= " << deg.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Degrader *>(deg.clone()));

}


gsell's avatar
gsell committed
668 669 670 671 672 673 674 675 676 677
/**
 *
 *
 * @param drift
 */
void ParallelCyclotronTracker::visitDrift(const Drift &drift) {
    *gmsg << "In drift L= " << drift.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Drift *>(drift.clone()));
}

678 679
/**
 *
680 681
 *
 *  @param
682 683 684 685 686
 */
void ParallelCyclotronTracker::visitFlexibleCollimator(const FlexibleCollimator &) {

}

687 688 689
void ParallelCyclotronTracker::visitOffset(const Offset & off) {
    if (opalRing_m == NULL)
        throw OpalException(
690 691
                            "ParallelCylcotronTracker::visitOffset",
                            "Attempt to place an offset when Ring not defined");
692 693
    Offset* offNonConst = const_cast<Offset*>(&off);
    offNonConst->updateGeometry(opalRing_m->getNextPosition(),
694
                                opalRing_m->getNextNormal());
695 696 697
    opalRing_m->appendElement(off);
}

gsell's avatar
gsell committed
698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729
/**
 *
 *
 * @param marker
 */
void ParallelCyclotronTracker::visitMarker(const Marker &marker) {
    //   *gmsg << "In Marker; L= " << marker.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Marker *>(marker.clone()));
    // Do nothing.
}

/**
 *
 *
 * @param corr
 */
void ParallelCyclotronTracker::visitMonitor(const Monitor &corr) {
    //   *gmsg << "In Monitor; L= " << corr.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Monitor *>(corr.clone()));
    //   applyDrift(flip_s * corr.getElementLength());
}

/**
 *
 *
 * @param mult
 */
void ParallelCyclotronTracker::visitMultipole(const Multipole &mult) {
    *gmsg << "In Multipole; L= " << mult.getElementLength() << " however the element is missing " << endl;
    myElements.push_back(dynamic_cast<Multipole *>(mult.clone()));
}

ext-rogers_c's avatar
ext-rogers_c committed
730 731 732 733 734 735 736
/**
 *
 *
 * @param multT
 */
void ParallelCyclotronTracker::visitMultipoleT(const MultipoleT &multT) {
    *gmsg << "Adding MultipoleT" << endl;
737
    if (opalRing_m != NULL) {
ext-rogers_c's avatar
ext-rogers_c committed
738
        opalRing_m->appendElement(multT);
739
    } else {
ext-rogers_c's avatar
ext-rogers_c committed
740 741
        throw OpalException("ParallelCyclotronTracker::visitMultipoleT",
                            "Need to define a RINGDEFINITION to use MultipoleT element");
742
    }
ext-rogers_c's avatar
ext-rogers_c committed
743 744 745
    myElements.push_back(dynamic_cast<MultipoleT *>(multT.clone()));
}

746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793
/**
 *
 *
 * @param multTstraight
 */
void ParallelCyclotronTracker::visitMultipoleTStraight(const MultipoleTStraight &multTstraight) {
    *gmsg << "Adding MultipoleTStraight" << endl;
    if (opalRing_m != NULL) {
        opalRing_m->appendElement(multTstraight);
    } else {
        throw OpalException("ParallelCyclotronTracker::visitMultipoleTStraight",
                            "Need to define a RINGDEFINITION to use MultipoleTStraight element");
    }
    myElements.push_back(dynamic_cast<MultipoleTStraight *>(multTstraight.clone()));
}

/**
 *
 *
 * @param multTccurv
 */
void ParallelCyclotronTracker::visitMultipoleTCurvedConstRadius(const MultipoleTCurvedConstRadius &multTccurv) {
    *gmsg << "Adding MultipoleTCurvedConstRadius" << endl;
    if (opalRing_m != NULL) {
        opalRing_m->appendElement(multTccurv);
    } else {
        throw OpalException("ParallelCyclotronTracker::visitMultipoleTCurvedConstRadius",
                            "Need to define a RINGDEFINITION to use MultipoleTCurvedConstRadius element");
    }
    myElements.push_back(dynamic_cast<MultipoleTCurvedConstRadius *>(multTccurv.clone()));
}

/**
 *
 *
 * @param multTvcurv
 */
void ParallelCyclotronTracker::visitMultipoleTCurvedVarRadius(const MultipoleTCurvedVarRadius &multTvcurv) {
    *gmsg << "Adding MultipoleTCurvedVarRadius" << endl;
    if (opalRing_m != NULL) {
        opalRing_m->appendElement(multTvcurv);
    } else {
        throw OpalException("ParallelCyclotronTracker::visitMultipoleTCurvedVarRadius",
                            "Need to define a RINGDEFINITION to use MultipoleTCurvedVarRadius element");
    }
    myElements.push_back(dynamic_cast<MultipoleTCurvedVarRadius *>(multTvcurv.clone()));
}

gsell's avatar
gsell committed
794 795 796 797 798 799
/**
 *
 *
 * @param prob
 */
void ParallelCyclotronTracker::visitProbe(const Probe &prob) {
800
    *gmsg << "* ------------------------------  Probe ------------------------------" << endl;
801 802
    Probe *elptr = dynamic_cast<Probe *>(prob.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
803

804 805
    *gmsg << "* Name    = " << elptr->getName() << endl;

806
    double xstart = elptr->getXStart();
807
    *gmsg << "* XStart  = " << xstart << " [m]" << endl;
gsell's avatar
gsell committed
808

809
    double xend = elptr->getXEnd();
810
    *gmsg << "* XEnd    = " << xend << " [m]" << endl;
gsell's avatar
gsell committed
811

812
    double ystart = elptr->getYStart();
813
    *gmsg << "* YStart  = " << ystart << " [m]" << endl;
gsell's avatar
gsell committed
814

815
    double yend = elptr->getYEnd();
816
    *gmsg << "* YEnd    = " << yend << " [m]" << endl;
gsell's avatar
gsell committed
817 818

    // initialise, do nothing
frey_m's avatar
frey_m committed
819
    elptr->initialise(itsBunch_m);
gsell's avatar
gsell committed
820

snuverink_j's avatar
snuverink_j committed
821
    double BcParameter[8] = {};
822

823 824 825 826 827
    BcParameter[0] = xstart;
    BcParameter[1] = xend;
    BcParameter[2] = ystart;
    BcParameter[3] = yend;
    BcParameter[4] = 1 ; // width
gsell's avatar
gsell committed
828 829

    // store probe parameters in the list
830
    buildupFieldList(BcParameter, ElementBase::PROBE, elptr);
gsell's avatar
gsell committed
831 832 833 834 835 836 837 838 839 840 841 842
}

/**
 *
 *
 * @param bend
 */
void ParallelCyclotronTracker::visitRBend(const RBend &bend) {
    *gmsg << "In RBend; L= " << bend.getElementLength() << " however the element is missing " << endl;
    myElements.push_back(dynamic_cast<RBend *>(bend.clone()));
}

843 844 845 846 847 848
void ParallelCyclotronTracker::visitSBend3D(const SBend3D &bend) {
    *gmsg << "Adding SBend3D" << endl;
    if (opalRing_m != NULL)
        opalRing_m->appendElement(bend);
    else
        throw OpalException("ParallelCyclotronTracker::visitSBend3D",
849
                            "Need to define a RINGDEFINITION to use SBend3D element");
850 851
}

ext-rogers_c's avatar
ext-rogers_c committed
852 853
void ParallelCyclotronTracker::visitScalingFFAMagnet(const ScalingFFAMagnet &bend) {
    *gmsg << "Adding ScalingFFAMagnet" << endl;
ext-rogers_c's avatar
ext-rogers_c committed
854
    if (opalRing_m != NULL) {
855
        opalRing_m->appendElement(bend);
ext-rogers_c's avatar
ext-rogers_c committed
856
    } else {
ext-rogers_c's avatar
ext-rogers_c committed
857 858
        throw OpalException("ParallelCyclotronTracker::visitScalingFFAMagnet",
                            "Need to define a RINGDEFINITION to use ScalingFFAMagnet element");
ext-rogers_c's avatar
ext-rogers_c committed
859
    }
860 861
}

862 863 864 865 866 867
void ParallelCyclotronTracker::visitVariableRFCavity(const VariableRFCavity &cav) {
    *gmsg << "Adding Variable RF Cavity" << endl;
    if (opalRing_m != NULL)
        opalRing_m->appendElement(cav);
    else
        throw OpalException("ParallelCyclotronTracker::visitVariableRFCavity",
868
                            "Need to define a RINGDEFINITION to use VariableRFCavity element");
869 870 871 872 873 874 875 876 877 878
}

void ParallelCyclotronTracker::visitVariableRFCavityFringeField
                                  (const VariableRFCavityFringeField &cav) {
    *gmsg << "Adding Variable RF Cavity with Fringe Field" << endl;
    if (opalRing_m != NULL)
        opalRing_m->appendElement(cav);
    else
        throw OpalException("ParallelCyclotronTracker::visitVariableRFCavityFringeField",
                            "Need to define a RINGDEFINITION to use VariableRFCavity element");
879 880
}

881 882 883 884 885 886 887 888 889
void ParallelCyclotronTracker::visitVerticalFFAMagnet(const VerticalFFAMagnet &mag) {
    *gmsg << "Adding Vertical FFA Magnet" << endl;
    if (opalRing_m != NULL)
        opalRing_m->appendElement(mag);
    else
        throw OpalException("ParallelCyclotronTracker::visitVerticalFFAMagnet",
                            "Need to define a RINGDEFINITION to use VerticalFFAMagnet element");
}

gsell's avatar
gsell committed
890 891 892 893 894 895 896
/**
 *
 *
 * @param as
 */
void ParallelCyclotronTracker::visitRFCavity(const RFCavity &as) {

897
    *gmsg << "* ------------------------------ RFCavity ------------------------------" << endl;
898

899 900
    RFCavity *elptr = dynamic_cast<RFCavity *>(as.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
901

902
    if ( elptr->getComponentType() != "SINGLEGAP" ) {
gsell's avatar
gsell committed
903 904 905 906 907
        *gmsg << (elptr->getComponentType()) << endl;
        throw OpalException("ParallelCyclotronTracker::visitRFCavity",
                            "The ParallelCyclotronTracker can only play with cyclotron type RF system currently ...");
    }

908
    double rmin = elptr->getRmin();
gsell's avatar
gsell committed
909 910
    *gmsg << "* Minimal radius of cavity= " << rmin << " [mm]" << endl;

911
    double rmax = elptr->getRmax();
gsell's avatar
gsell committed
912 913
    *gmsg << "* Maximal radius of cavity= " << rmax << " [mm]" << endl;

914
    double rff = elptr->getCycFrequency();
gsell's avatar
gsell committed
915 916
    *gmsg << "* RF frequency (2*pi*f)= " << rff << " [rad/s]" << endl;

917
    std::string fmfn = elptr->getFieldMapFN();
gsell's avatar
gsell committed
918
    *gmsg << "* RF Field map file name= " << fmfn << endl;
919

920
    double angle = elptr->getAzimuth();
gsell's avatar
gsell committed
921 922
    *gmsg << "* Cavity azimuth position= " << angle << " [deg] " << endl;

923
    double gap = elptr->getGapWidth();
gsell's avatar
gsell committed
924 925
    *gmsg << "* Cavity gap width= " << gap << " [mm] " << endl;

926
    double pdis = elptr->getPerpenDistance();
gsell's avatar
gsell committed
927 928
    *gmsg << "* Cavity Shift distance= " << pdis << " [mm] " << endl;

929
    double phi0 = elptr->getPhi0();
gsell's avatar
gsell committed
930 931
    *gmsg << "* Initial RF phase (t=0)= " << phi0 << " [deg] " << endl;

932 933 934
    /*
      Setup time dependence and in case of no
      timedependence use a polynom with  a_0 = 1 and a_k = 0, k = 1,2,3.
935
    */
936

937 938 939 940
    std::shared_ptr<AbstractTimeDependence> freq_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> ampl_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> phase_atd = nullptr;

941
    dvector_t  unityVec;
942 943 944 945
    unityVec.push_back(1.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
946

947
    if (elptr->getFrequencyModelName() != "") {
948 949
        freq_atd = AbstractTimeDependence::getTimeDependence(elptr->getFrequencyModelName());
        *gmsg << "* Variable frequency RF Model name " << elptr->getFrequencyModelName() << endl;
950
    }
951 952
    else
        freq_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
953 954

    if (elptr->getAmplitudeModelName() != "") {
955 956
        ampl_atd = AbstractTimeDependence::getTimeDependence(elptr->getAmplitudeModelName());
        *gmsg << "* Variable amplitude RF Model name " << elptr->getAmplitudeModelName() << endl;
957
    }
958 959
    else
        ampl_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
960 961

    if (elptr->getPhaseModelName() != "") {
962 963
        phase_atd = AbstractTimeDependence::getTimeDependence(elptr->getPhaseModelName());
        *gmsg << "* Variable phase RF Model name " << elptr->getPhaseModelName() << endl;
964
    }
965 966
    else
        phase_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
967

gsell's avatar
gsell committed
968
    // read cavity voltage profile data from file.
frey_m's avatar
frey_m committed
969
    elptr->initialise(itsBunch_m, freq_atd, ampl_atd, phase_atd);
gsell's avatar
gsell committed
970

snuverink_j's avatar
snuverink_j committed
971
    double BcParameter[8] = {};
972

973 974
    BcParameter[0] = 0.001 * rmin;
    BcParameter[1] = 0.001 * rmax;
975
    BcParameter[2] = 0.001 * pdis;
gsell's avatar
gsell committed
976 977
    BcParameter[3] = angle;

978
    buildupFieldList(BcParameter, ElementBase::RFCAVITY, elptr);
gsell's avatar
gsell committed
979 980 981 982 983 984 985 986
}

/**
 *
 *
 * @param bend
 */
void ParallelCyclotronTracker::visitSBend(const SBend &bend) {
987
    *gmsg << "In SBend; L = " << bend.getElementLength() << " however the element is missing " << endl;
gsell's avatar
gsell committed
988 989 990 991 992 993 994 995 996
    myElements.push_back(dynamic_cast<SBend *>(bend.clone()));
}

/**
 *
 *
 * @param sept
 */
void ParallelCyclotronTracker::visitSeptum(const Septum &sept) {
997

998
    *gmsg << endl << "* ------------------------------ Septum ------------------------------- *" << endl;
gsell's avatar
gsell committed
999

1000 1001
    Septum *elptr = dynamic_cast<Septum *>(sept.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
1002

1003
    double xstart = elptr->getXStart();
1004
    *gmsg << "* XStart  = " << xstart << " [m]" << endl;
gsell's avatar
gsell committed
1005

1006
    double xend = elptr->getXEnd();
1007
    *gmsg << "* XEnd    = " << xend << " [m]" << endl;
gsell's avatar
gsell committed
1008

1009
    double ystart = elptr->getYStart();
1010
    *gmsg << "* YStart  = " << ystart << " [m]" << endl;
gsell's avatar
gsell committed
1011

1012
    double yend = elptr->getYEnd();
1013
    *gmsg << "* YEnd    = " << yend << " [m]" << endl;
gsell's avatar
gsell committed
1014

1015
    double width = elptr->getWidth();
1016
    *gmsg << "* Width   = " << width << " [m]" << endl;
gsell's avatar
gsell committed
1017 1018 1019


    // initialise, do nothing
frey_m's avatar
frey_m committed
1020
    elptr->initialise(itsBunch_m);
gsell's avatar
gsell committed
1021

snuverink_j's avatar
snuverink_j committed
1022
    double BcParameter[8] = {};
1023

1024 1025 1026 1027 1028
    BcParameter[0] = xstart;
    BcParameter[1] = xend;
    BcParameter[2] = ystart;
    BcParameter[3] = yend;
    BcParameter[4] = width;
gsell's avatar
gsell committed
1029 1030

    // store septum parameters in the list
1031
    buildupFieldList(BcParameter, ElementBase::SEPTUM, elptr);
gsell's avatar
gsell committed
1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047
}

/**
 *
 *
 * @param solenoid
 */
void ParallelCyclotronTracker::visitSolenoid(const Solenoid &solenoid) {
    myElements.push_back(dynamic_cast<Solenoid *>(solenoid.clone()));
    Component *elptr = *(--myElements.end());
    if(!elptr->hasAttribute("ELEMEDGE")) {
        *gmsg << "Solenoid: no position of the element given!" << endl;
        return;
    }
}

1048 1049 1050 1051 1052
/**
 *
 *
 * @param stripper
 */
1053

1054
void ParallelCyclotronTracker::visitStripper(const Stripper &stripper) {
1055

1056
    *gmsg << "* ------------------------------ Stripper ------------------------------" << endl;
1057