ParallelCyclotronTracker.cpp 123 KB
Newer Older
gsell's avatar
gsell committed
1 2 3 4 5 6 7 8 9
// ------------------------------------------------------------------------
// $RCSfile: ParallelCyclotronTracker.cpp,v $
// ------------------------------------------------------------------------
// $Revision: 1.1 $initialLocalNum_m
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: ParallelCyclotronTracker
ext-rogers_c's avatar
ext-rogers_c committed
10
//   The class for tracking particles with 3D space charge in Cyclotrons and FFAs
gsell's avatar
gsell committed
11 12 13 14
//
// ------------------------------------------------------------------------
//
// $Date: 2007/10/17 04:00:08 $
15
// $Author: adelmann, yang, winklehner $
gsell's avatar
gsell committed
16 17
//
// ------------------------------------------------------------------------
kraus's avatar
kraus committed
18 19

#include "Algorithms/ParallelCyclotronTracker.h"
20

gsell's avatar
gsell committed
21
#include <fstream>
22 23
#include <iostream>
#include <limits>
gsell's avatar
gsell committed
24
#include <vector>
frey_m's avatar
frey_m committed
25
#include <numeric>
26 27

#include "AbstractObjects/Element.h"
28
#include "AbstractObjects/OpalData.h"
gsell's avatar
gsell committed
29

30
#include "AbsBeamline/CCollimator.h"
gsell's avatar
gsell committed
31 32
#include "AbsBeamline/Corrector.h"
#include "AbsBeamline/Cyclotron.h"
33
#include "AbsBeamline/CyclotronValley.h"
adelmann's avatar
adelmann committed
34
#include "AbsBeamline/Degrader.h"
gsell's avatar
gsell committed
35 36 37 38
#include "AbsBeamline/Diagnostic.h"
#include "AbsBeamline/Drift.h"
#include "AbsBeamline/ElementBase.h"
#include "AbsBeamline/Lambertson.h"
39
#include "AbsBeamline/Offset.h"
gsell's avatar
gsell committed
40 41 42
#include "AbsBeamline/Marker.h"
#include "AbsBeamline/Monitor.h"
#include "AbsBeamline/Multipole.h"
ext-rogers_c's avatar
ext-rogers_c committed
43
#include "AbsBeamline/MultipoleT.h"
44 45 46 47
#include "AbsBeamline/MultipoleTBase.h"
#include "AbsBeamline/MultipoleTStraight.h"
#include "AbsBeamline/MultipoleTCurvedConstRadius.h"
#include "AbsBeamline/MultipoleTCurvedVarRadius.h"
48
#include "AbsBeamline/PluginElement.h"
gsell's avatar
gsell committed
49 50
#include "AbsBeamline/Probe.h"
#include "AbsBeamline/RBend.h"
51
#include "AbsBeamline/Ring.h"
gsell's avatar
gsell committed
52 53 54
#include "AbsBeamline/RFCavity.h"
#include "AbsBeamline/RFQuadrupole.h"
#include "AbsBeamline/SBend.h"
55
#include "AbsBeamline/SBend3D.h"
ext-rogers_c's avatar
ext-rogers_c committed
56
#include "AbsBeamline/ScalingFFAMagnet.h"
gsell's avatar
gsell committed
57 58 59 60
#include "AbsBeamline/Separator.h"
#include "AbsBeamline/Septum.h"
#include "AbsBeamline/Solenoid.h"
#include "AbsBeamline/Stripper.h"
61
#include "AbsBeamline/VariableRFCavity.h"
62
#include "AbsBeamline/VariableRFCavityFringeField.h"
63

64 65
#include "Algorithms/Ctunes.h"
#include "Algorithms/PolynomialTimeDependence.h"
gsell's avatar
gsell committed
66 67

#include "Beamlines/Beamline.h"
68
#include "Beamlines/FlaggedBeamline.h"
gsell's avatar
gsell committed
69

70
#include "Elements/OpalBeamline.h"
gsell's avatar
gsell committed
71 72 73 74

#include "Physics/Physics.h"

#include "Utilities/OpalException.h"
75
#include "Utilities/Options.h"
gsell's avatar
gsell committed
76

77
#include "BasicActions/DumpFields.h"
78
#include "BasicActions/DumpEMFields.h"
79

80
#include "Structure/BoundaryGeometry.h"
81
#include "Structure/DataSink.h"
frey_m's avatar
frey_m committed
82

gsell's avatar
gsell committed
83 84 85
using Physics::pi;
using Physics::q_e;

86

87
const double c_mmtns = Physics::c * 1.0e-6; // m/s --> mm/ns
88

89 90 91
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
92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107

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
108
                                                   PartBunchBase<double, 3> *bunch,
gsell's avatar
gsell committed
109 110 111
                                                   DataSink &ds,
                                                   const PartData &reference,
                                                   bool revBeam, bool revTrack,
frey_m's avatar
frey_m committed
112
                                                   int maxSTEPS, int timeIntegrator,
113 114 115 116 117
                                                   const int& numBunch,
                                                   const double& mbEta,
                                                   const double& mbPara,
                                                   const std::string& mbMode,
                                                   const std::string& mbBinning)
frey_m's avatar
frey_m committed
118 119 120 121 122 123 124 125 126 127
    : 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
128 129 130
    itsBeamline = dynamic_cast<Beamline *>(beamline.clone());
    itsDataSink = &ds;

131 132 133
    if ( numBunch > 1 ) {
        mbHandler_m = std::unique_ptr<MultiBunchHandler>(
            new MultiBunchHandler(bunch, numBunch, mbEta,
frey_m's avatar
frey_m committed
134
                                  mbPara, mbMode, mbBinning)
135 136 137
        );
    }

gsell's avatar
gsell committed
138 139 140 141
    IntegrationTimer_m = IpplTimings::getTimer("Integration");
    TransformTimer_m   = IpplTimings::getTimer("Frametransform");
    DumpTimer_m        = IpplTimings::getTimer("Dump");
    BinRepartTimer_m   = IpplTimings::getTimer("Binaryrepart");
142 143
    PluginElemTimer_m  = IpplTimings::getTimer("PluginElements");
    DelParticleTimer_m = IpplTimings::getTimer("DeleteParticles");
144

145 146 147 148 149 150 151 152 153
    // 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;
154

155 156 157 158 159 160 161 162
    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
163 164 165 166 167 168 169
}

/**
 * Destructor ParallelCyclotronTracker
 *
 */
ParallelCyclotronTracker::~ParallelCyclotronTracker() {
170 171
    for(Component* component : myElements) {
        delete(component);
gsell's avatar
gsell committed
172
    }
173 174
    for(auto fd : FieldDimensions) {
        delete(fd);
gsell's avatar
gsell committed
175 176
    }
    delete itsBeamline;
177
    // delete opalRing_m;
gsell's avatar
gsell committed
178 179
}

frey_m's avatar
frey_m committed
180

181 182 183 184 185 186
/**
 *
 *
 * @param fn Base file name
 */
void ParallelCyclotronTracker::bgf_main_collision_test() {
187
    if(!bgf_m) return;
188

189
    Inform msg("bgf_main_collision_test ");
190

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

195
    Vector_t intecoords = 0.0;
196

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

200
    int triId = 0;
frey_m's avatar
frey_m committed
201 202 203
    for(size_t i = 0; i < itsBunch_m->getLocalNum(); i++) {
        int res = bgf_m->partInside(itsBunch_m->R[i], itsBunch_m->P[i], dtime, itsBunch_m->PType[i], itsBunch_m->Q[i], intecoords, triId);
        //int res = bgf_m->partInside(itsBunch_m->R[i]*1.0e-3, itsBunch_m->P[i], dtime, itsBunch_m->PType[i], itsBunch_m->Q[i], intecoords, triId);
204
        if(res >= 0) {
frey_m's avatar
frey_m committed
205
            itsBunch_m->Bin[i] = -1;
206
        }
207
    }
208 209
}

210
// only used for dumping into stat file
frey_m's avatar
frey_m committed
211 212 213 214 215
void ParallelCyclotronTracker::dumpAngle(const double& theta,
                                         double& prevAzimuth,
                                         double& azimuth)
{
    if ( prevAzimuth < 0.0 ) { // only at first occurrence
frey_m's avatar
frey_m committed
216 217 218 219 220
        double plus = 0.0;
        if ( OpalData::getInstance()->inRestartRun() ) {
            plus = 360.0 * turnnumber_m;
        }
        azimuth_m = theta + plus;
221
    } else {
frey_m's avatar
frey_m committed
222
        double dtheta = theta - prevAzimuth;
223 224 225 226 227 228
        if ( dtheta < 0 ) {
            dtheta += 360.0;
        }
        if ( dtheta > 180 ) { // rotating clockwise, reduce angle
            dtheta -= 360;
        }
frey_m's avatar
frey_m committed
229
        azimuth += dtheta;
230
    }
frey_m's avatar
frey_m committed
231
    prevAzimuth = theta;
232 233 234
}


frey_m's avatar
frey_m committed
235 236 237 238 239 240
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
241 242
void ParallelCyclotronTracker::computePathLengthUpdate(std::vector<double>& dl,
                                                       const double& dt)
243
{
frey_m's avatar
frey_m committed
244 245
    // the last element in dotP is the dot-product over all particles
    std::vector<double> dotP(dl.size());
246
    if ( Options::psDumpFrame == Options::BUNCH_MEAN || isMultiBunch()) {
247

frey_m's avatar
frey_m committed
248 249
        for(unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) {
            dotP[itsBunch_m->bunchNum[i]] += dot(itsBunch_m->P[i], itsBunch_m->P[i]);
250 251
        }

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

frey_m's avatar
frey_m committed
254 255 256
        // 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
257

frey_m's avatar
frey_m committed
258 259 260 261
        // bunch specific --> multi-bunches only
        for (short b = 0; b < (short)dotP.size() - 1; ++b) {
            dotP[b] /= double(itsBunch_m->getTotalNumPerBunch(b));
        }
262 263

    } else if ( itsBunch_m->getLocalNum() == 0 ) {
264
        // here we are in Options::GLOBAL mode
frey_m's avatar
frey_m committed
265
        dotP[0] = 0.0;
266
    } else {
267
        // here we are in Options::GLOBAL mode
frey_m's avatar
frey_m committed
268
        dotP[0] = dot(itsBunch_m->P[0], itsBunch_m->P[0]);
269 270
    }

frey_m's avatar
frey_m committed
271 272 273 274 275
    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
    }
276 277 278
}


279

gsell's avatar
gsell committed
280 281 282 283 284
/**
 *
 *
 * @param fn Base file name
 */
285
void ParallelCyclotronTracker::openFiles(size_t numFiles, std::string SfileName) {
gsell's avatar
gsell committed
286

287
    for (unsigned int i=0; i<numFiles; i++) {
288 289 290 291 292 293 294 295
        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
296

297 298 299 300 301 302
        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;
303
    }
gsell's avatar
gsell committed
304 305 306 307 308 309 310 311
}

/**
 * Close all files related to
 * special output in the Cyclotron
 * mode.
 */
void ParallelCyclotronTracker::closeFiles() {
312
    for (auto & file : outfTheta_m) {
313
        file->close();
314
    }
gsell's avatar
gsell committed
315 316
}

317
/**
318 319 320
 *
 * @param ring
 */
kraus's avatar
kraus committed
321
void ParallelCyclotronTracker::visitRing(const Ring &ring) {
Daniel Winklehner's avatar
Daniel Winklehner committed
322

323
    *gmsg << "* ----------------------------- Adding Ring ------------------------------ *" << endl;
Daniel Winklehner's avatar
Daniel Winklehner committed
324

snuverink_j's avatar
snuverink_j committed
325
    delete opalRing_m;
Daniel Winklehner's avatar
Daniel Winklehner committed
326

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

329
    myElements.push_back(opalRing_m);
Daniel Winklehner's avatar
Daniel Winklehner committed
330

frey_m's avatar
frey_m committed
331
    opalRing_m->initialise(itsBunch_m);
332 333 334 335

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

337
    if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
kraus's avatar
kraus committed
338
        throw OpalException("Error in ParallelCyclotronTracker::visitRing",
339 340
                            "PHIINIT is out of [-180, 180)!");
    }
Daniel Winklehner's avatar
Daniel Winklehner committed
341 342

    referenceZ = 0.0;
343
    referencePz = 0.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
344 345 346 347

    referencePtot = itsReference.getGamma() * itsReference.getBeta();
    referencePt = sqrt(referencePtot * referencePtot - referencePr * referencePr);

348 349
    if(referencePtot < 0.0)
        referencePt *= -1.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
350

351 352
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
gsell's avatar
gsell committed
353

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

356
    buildupFieldList(BcParameter, ElementBase::RING, opalRing_m);
357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372

    // 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
373 374 375 376 377 378 379 380

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

381
    *gmsg << "* -------------------------- Adding Cyclotron ---------------------------- *" << endl;
gsell's avatar
gsell committed
382

383 384
    Cyclotron *elptr = dynamic_cast<Cyclotron *>(cycl.clone());
    myElements.push_back(elptr);
385

386
    // Is this a Spiral Inflector Simulation? If yes, we'll give the user some
387 388 389 390 391 392
    // useful information
    spiral_flag = elptr->getSpiralFlag();

    if(spiral_flag) {

        *gmsg << endl << "* This is a Spiral Inflector Simulation! This means the following:" << endl;
393
        *gmsg         << "* 1.) It is up to the user to provide appropriate geometry, electric and magnetic fields!" << endl;
394 395
        *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
396
        *gmsg         << "* 2.) For high currents it is strongly recommended to use the SAAMG fieldsolver," << endl;
397 398
        *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;
399 400
        *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;
401 402 403
        *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
404 405 406
        if (isMultiBunch()) {
            mbHandler_m = nullptr;
        }
407 408
    }

409
    // Fresh run (no restart):
410
    if(!OpalData::getInstance()->inRestartRun()) {
411

412
        // Get reference values from cyclotron element
413
        // For now, these are still stored in mm. should be the only ones. -DW
414 415 416 417
        referenceR     = elptr->getRinit();
        referenceTheta = elptr->getPHIinit();
        referenceZ     = elptr->getZinit();
        referencePr    = elptr->getPRinit();
418 419
        referencePz    = elptr->getPZinit();

420
        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
421 422
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                "PHIINIT is out of [-180, 180)!");
423 424 425 426 427
        }

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

        // Calculate reference azimuthal (tangential) momentum from total-, z- and radial momentum:
428
        float insqrt = referencePtot * referencePtot - \
429
            referencePr * referencePr - referencePz * referencePz;
430 431 432 433 434

        if(insqrt < 0) {

            if(insqrt > -1.0e-10) {

435
                referencePt = 0.0;
436 437 438

            } else {

439
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
440
                                    "Pt imaginary!");
441 442 443 444 445 446 447
            }

        } else {

            referencePt = sqrt(insqrt);
        }

448
        if(referencePtot < 0.0)
Daniel Winklehner's avatar
Daniel Winklehner committed
449
            referencePt *= -1.0;
450 451
        // End calculate referencePt

452
        // Restart a run:
453 454
    } else {

455
        // If the user wants to save the restarted run in local frame,
456
        // make sure the previous h5 file was local too
457
      if (Options::psDumpFrame != Options::GLOBAL) {
winklehner_d's avatar
-DW  
winklehner_d committed
458
        if (!previousH5Local) {
459 460
                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
461
        }
462 463
            // 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
464 465
      } else {
        if (previousH5Local) {
466 467
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                    "You are trying a global restart from a local h5 file!");
468
            }
469
        }
470

471
        // Adjust some of the reference variables from the h5 file
472 473
        referencePhi *= Physics::deg2rad;
        referencePsi *= Physics::deg2rad;
474
        referencePtot = bega;
475
        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
476
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
477 478
                                "PHIINIT is out of [-180, 180)!");
        }
479 480
    }

481 482
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
483

484
    *gmsg << endl;
adelmann's avatar
adelmann committed
485
    *gmsg << "* Bunch global starting position:" << endl;
486 487
    *gmsg << "* RINIT = " << referenceR  << " [mm]" << endl;
    *gmsg << "* PHIINIT = " << referenceTheta << " [deg]" << endl;
488
    *gmsg << "* ZINIT = " << referenceZ << " [mm]" << endl;
489
    *gmsg << endl;
adelmann's avatar
adelmann committed
490
    *gmsg << "* Bunch global starting momenta:" << endl;
gsell's avatar
gsell committed
491 492
    *gmsg << "* Initial gamma = " << itsReference.getGamma() << endl;
    *gmsg << "* Initial beta = " << itsReference.getBeta() << endl;
493
    *gmsg << "* Reference total momentum (beta * gamma) = " << referencePtot * 1000.0 << " [MCU]" << endl;
494 495 496
    *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;
497
    *gmsg << endl;
adelmann's avatar
adelmann committed
498

gsell's avatar
gsell committed
499
    double sym = elptr->getSymmetry();
snuverink_j's avatar
snuverink_j committed
500
    *gmsg << "* " << sym << "-fold field symmetry " << endl;
gsell's avatar
gsell committed
501

502 503 504
    // ckr: this just returned the default value as defined in Component.h
    // double rff = elptr->getRfFrequ();
    // *gmsg << "* Rf frequency= " << rff << " [MHz]" << endl;
gsell's avatar
gsell committed
505

506
    std::string fmfn = elptr->getFieldMapFN();
507
    *gmsg << "* Field map file name = " << fmfn << " " << endl;
gsell's avatar
gsell committed
508

509
    std::string type = elptr->getCyclotronType();
510
    *gmsg << "* Type of cyclotron = " << type << " " << endl;
511

512 513
    double rmin = elptr->getMinR();
    double rmax = elptr->getMaxR();
adelmann's avatar
adelmann committed
514
    *gmsg << "* Radial aperture = " << rmin << " ... " << rmax<<" [m] "<< endl;
515 516 517

    double zmin = elptr->getMinZ();
    double zmax = elptr->getMaxZ();
adelmann's avatar
adelmann committed
518
    *gmsg << "* Vertical aperture = " << zmin << " ... " << zmax<<" [m]"<< endl;
gsell's avatar
gsell committed
519 520

    double h = elptr->getCyclHarm();
521
    *gmsg << "* Number of trimcoils = " << elptr->getNumberOfTrimcoils() << endl;
522
    *gmsg << "* Harmonic number h = " << h << " " << endl;
523

gsell's avatar
gsell committed
524 525 526 527 528 529 530
    /**
     * 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
531
     * fieldflag = 5, readin FFA format file for MSU/FNAL FFA
gsell's avatar
gsell committed
532
     * fieldflag = 6, readin both median plane B field map and 3D E field map of RF cavity for compact cyclotron
533
     * fieldflag = 7, read in fields for Daniel's synchrocyclotron simulations
gsell's avatar
gsell committed
534
     */
535
    int fieldflag = elptr->getFieldFlag(type);
gsell's avatar
gsell committed
536

537
    // Read in cyclotron field maps (midplane + 3D fields if desired).
frey_m's avatar
frey_m committed
538
    elptr->initialise(itsBunch_m, fieldflag, elptr->getBScale());
gsell's avatar
gsell committed
539

snuverink_j's avatar
snuverink_j committed
540
    double BcParameter[8] = {};
541

542 543
    BcParameter[0] = 0.001 * elptr->getRmin();
    BcParameter[1] = 0.001 * elptr->getRmax();
gsell's avatar
gsell committed
544

545
    // Store inner radius and outer radius of cyclotron field map in the list
546
    buildupFieldList(BcParameter, ElementBase::CYCLOTRON, elptr);
gsell's avatar
gsell committed
547 548 549 550 551 552 553 554 555 556 557 558 559 560 561
}

/**
 * Not implemented and most probable never used
 *
 */
void ParallelCyclotronTracker::visitBeamBeam(const BeamBeam &) {
    *gmsg << "In BeamBeam tracker is missing " << endl;
}

/**
 *
 *
 * @param coll
 */
562
void ParallelCyclotronTracker::visitCCollimator(const CCollimator &coll) {
gsell's avatar
gsell committed
563

564
    *gmsg << "* --------- Collimator -----------------------------" << endl;
gsell's avatar
gsell committed
565

566
    CCollimator* elptr = dynamic_cast<CCollimator *>(coll.clone());
567
    myElements.push_back(elptr);
gsell's avatar
gsell committed
568

569
    double xstart = elptr->getXStart();
adelmann's avatar
adelmann committed
570
    *gmsg << "* Xstart= " << xstart << " [mm]" << endl;
gsell's avatar
gsell committed
571

572
    double xend = elptr->getXEnd();
adelmann's avatar
adelmann committed
573
    *gmsg << "* Xend= " << xend << " [mm]" << endl;
gsell's avatar
gsell committed
574

575
    double ystart = elptr->getYStart();
adelmann's avatar
adelmann committed
576
    *gmsg << "* Ystart= " << ystart << " [mm]" << endl;
gsell's avatar
gsell committed
577

578
    double yend = elptr->getYEnd();
579
    *gmsg << "* Yend= " << yend << " [mm]" << endl;
gsell's avatar
gsell committed
580

581
    double zstart = elptr->getZStart();
adelmann's avatar
adelmann committed
582
    *gmsg << "* Zstart= " << zstart << " [mm]" << endl;
583 584

    double zend = elptr->getZEnd();
585
    *gmsg << "* Zend= " << zend << " [mm]" << endl;
586

587
    double width = elptr->getWidth();
adelmann's avatar
adelmann committed
588
    *gmsg << "* Width= " << width << " [mm]" << endl;
gsell's avatar
gsell committed
589

frey_m's avatar
frey_m committed
590
    elptr->initialise(itsBunch_m);
gsell's avatar
gsell committed
591

snuverink_j's avatar
snuverink_j committed
592
    double BcParameter[8] = {};
593

594 595 596 597 598 599
    BcParameter[0] = 0.001 * xstart ;
    BcParameter[1] = 0.001 * xend;
    BcParameter[2] = 0.001 * ystart ;
    BcParameter[3] = 0.001 * yend;
    BcParameter[4] = 0.001 * width ;

600
    buildupFieldList(BcParameter, ElementBase::CCOLLIMATOR, elptr);
gsell's avatar
gsell committed
601 602 603 604 605 606 607 608 609 610 611 612
}

/**
 *
 *
 * @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
613 614 615 616 617 618 619 620 621 622 623 624
/**
 *
 *
 * @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
625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644
/**
 *
 *
 * @param diag
 */
void ParallelCyclotronTracker::visitDiagnostic(const Diagnostic &diag) {
    *gmsg << "In Diagnostic; L= " << diag.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Diagnostic *>(diag.clone()));
}

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

645 646
/**
 *
647 648
 *
 *  @param
649 650 651 652 653
 */
void ParallelCyclotronTracker::visitFlexibleCollimator(const FlexibleCollimator &) {

}

gsell's avatar
gsell committed
654 655 656 657 658 659 660 661 662 663
/**
 *
 *
 * @param lamb
 */
void ParallelCyclotronTracker::visitLambertson(const Lambertson &lamb) {
    *gmsg << "In Lambertson; L= " << lamb.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Lambertson *>(lamb.clone()));
}

664 665 666
void ParallelCyclotronTracker::visitOffset(const Offset & off) {
    if (opalRing_m == NULL)
        throw OpalException(
667 668
                            "ParallelCylcotronTracker::visitOffset",
                            "Attempt to place an offset when Ring not defined");
669 670
    Offset* offNonConst = const_cast<Offset*>(&off);
    offNonConst->updateGeometry(opalRing_m->getNextPosition(),
671
                                opalRing_m->getNextNormal());
672 673 674
    opalRing_m->appendElement(off);
}

gsell's avatar
gsell committed
675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706
/**
 *
 *
 * @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
707 708 709 710 711 712 713
/**
 *
 *
 * @param multT
 */
void ParallelCyclotronTracker::visitMultipoleT(const MultipoleT &multT) {
    *gmsg << "Adding MultipoleT" << endl;
714
    if (opalRing_m != NULL) {
ext-rogers_c's avatar
ext-rogers_c committed
715
        opalRing_m->appendElement(multT);
716
    } else {
ext-rogers_c's avatar
ext-rogers_c committed
717 718
        throw OpalException("ParallelCyclotronTracker::visitMultipoleT",
                            "Need to define a RINGDEFINITION to use MultipoleT element");
719
    }
ext-rogers_c's avatar
ext-rogers_c committed
720 721 722
    myElements.push_back(dynamic_cast<MultipoleT *>(multT.clone()));
}

723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 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
/**
 *
 *
 * @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
771 772 773 774 775 776
/**
 *
 *
 * @param prob
 */
void ParallelCyclotronTracker::visitProbe(const Probe &prob) {
777
    *gmsg << "* -----------  Probe -------------------------------" << endl;
778 779
    Probe *elptr = dynamic_cast<Probe *>(prob.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
780

781
    double xstart = elptr->getXStart();
782
    *gmsg << "XStart= " << xstart << " [mm]" << endl;
gsell's avatar
gsell committed
783

784
    double xend = elptr->getXEnd();
785
    *gmsg << "XEnd= " << xend << " [mm]" << endl;
gsell's avatar
gsell committed
786

787
    double ystart = elptr->getYStart();
788
    *gmsg << "YStart= " << ystart << " [mm]" << endl;
gsell's avatar
gsell committed
789

790
    double yend = elptr->getYEnd();
791
    *gmsg << "YEnd= " << yend << " [mm]" << endl;
gsell's avatar
gsell committed
792 793

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

snuverink_j's avatar
snuverink_j committed
796
    double BcParameter[8] = {};
797

798 799 800 801
    BcParameter[0] = 0.001 * xstart ;
    BcParameter[1] = 0.001 * xend;
    BcParameter[2] = 0.001 * ystart ;
    BcParameter[3] = 0.001 * yend;
802
    BcParameter[4] = 0.001 ; // width
gsell's avatar
gsell committed
803 804

    // store probe parameters in the list
805
    buildupFieldList(BcParameter, ElementBase::PROBE, elptr);
gsell's avatar
gsell committed
806 807 808 809 810 811 812 813 814 815 816 817
}

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

818 819 820 821 822 823
void ParallelCyclotronTracker::visitSBend3D(const SBend3D &bend) {
    *gmsg << "Adding SBend3D" << endl;
    if (opalRing_m != NULL)
        opalRing_m->appendElement(bend);
    else
        throw OpalException("ParallelCyclotronTracker::visitSBend3D",
824
                            "Need to define a RINGDEFINITION to use SBend3D element");
825 826
}

ext-rogers_c's avatar
ext-rogers_c committed
827 828
void ParallelCyclotronTracker::visitScalingFFAMagnet(const ScalingFFAMagnet &bend) {
    *gmsg << "Adding ScalingFFAMagnet" << endl;
ext-rogers_c's avatar
ext-rogers_c committed
829
    if (opalRing_m != NULL) {
830
        opalRing_m->appendElement(bend);
ext-rogers_c's avatar
ext-rogers_c committed
831
    } else {
ext-rogers_c's avatar
ext-rogers_c committed
832 833
        throw OpalException("ParallelCyclotronTracker::visitScalingFFAMagnet",
                            "Need to define a RINGDEFINITION to use ScalingFFAMagnet element");
ext-rogers_c's avatar
ext-rogers_c committed
834
    }
835 836
}

837 838 839 840 841 842
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",
843
                            "Need to define a RINGDEFINITION to use VariableRFCavity element");
844 845 846 847 848 849 850 851 852 853
}

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");
854 855
}

gsell's avatar
gsell committed
856 857 858 859 860 861 862 863
/**
 *
 *
 * @param as
 */
void ParallelCyclotronTracker::visitRFCavity(const RFCavity &as) {

    *gmsg << "* --------- RFCavity ------------------------------" << endl;
864

865 866
    RFCavity *elptr = dynamic_cast<RFCavity *>(as.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
867 868 869 870 871 872 873

    if((elptr->getComponentType() != "SINGLEGAP") && (elptr->getComponentType() != "DOUBLEGAP")) {
        *gmsg << (elptr->getComponentType()) << endl;
        throw OpalException("ParallelCyclotronTracker::visitRFCavity",
                            "The ParallelCyclotronTracker can only play with cyclotron type RF system currently ...");
    }

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

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

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

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

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

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

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

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

898 899 900
    /*
      Setup time dependence and in case of no
      timedependence use a polynom with  a_0 = 1 and a_k = 0, k = 1,2,3.
901
    */
902

903 904 905 906
    std::shared_ptr<AbstractTimeDependence> freq_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> ampl_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> phase_atd = nullptr;

907
    dvector_t  unityVec;
908 909 910 911
    unityVec.push_back(1.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
912

913
    if (elptr->getFrequencyModelName() != "") {