ParallelCyclotronTracker.cpp 126 KB
Newer Older
gsell's avatar
gsell committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14
// ------------------------------------------------------------------------
// $RCSfile: ParallelCyclotronTracker.cpp,v $
// ------------------------------------------------------------------------
// $Revision: 1.1 $initialLocalNum_m
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: ParallelCyclotronTracker
//   The class for tracking particles with 3D space charge in Cyclotrons and FFAG's
//
// ------------------------------------------------------------------------
//
// $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/Ctunes.hpp>
kraus's avatar
kraus committed
20
#include "Algorithms/ParallelCyclotronTracker.h"
21 22
#include "Algorithms/PolynomialTimeDependence.h"
#include "Elements/OpalPolynomialTimeDependence.h"
23

gsell's avatar
gsell committed
24 25 26 27
#include <cfloat>
#include <iostream>
#include <fstream>
#include <vector>
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"
adelmann's avatar
adelmann committed
33
#include "AbsBeamline/Degrader.h"
gsell's avatar
gsell committed
34 35 36 37
#include "AbsBeamline/Diagnostic.h"
#include "AbsBeamline/Drift.h"
#include "AbsBeamline/ElementBase.h"
#include "AbsBeamline/Lambertson.h"
38
#include "AbsBeamline/Offset.h"
gsell's avatar
gsell committed
39 40 41
#include "AbsBeamline/Marker.h"
#include "AbsBeamline/Monitor.h"
#include "AbsBeamline/Multipole.h"
ext-rogers_c's avatar
ext-rogers_c committed
42
#include "AbsBeamline/MultipoleT.h"
43 44 45 46
#include "AbsBeamline/MultipoleTBase.h"
#include "AbsBeamline/MultipoleTStraight.h"
#include "AbsBeamline/MultipoleTCurvedConstRadius.h"
#include "AbsBeamline/MultipoleTCurvedVarRadius.h"
gsell's avatar
gsell committed
47 48 49 50 51
#include "AbsBeamline/Probe.h"
#include "AbsBeamline/RBend.h"
#include "AbsBeamline/RFCavity.h"
#include "AbsBeamline/RFQuadrupole.h"
#include "AbsBeamline/SBend.h"
52
#include "AbsBeamline/SBend3D.h"
53
#include "AbsBeamline/ScalingFFAGMagnet.h"
gsell's avatar
gsell committed
54 55 56 57 58
#include "AbsBeamline/Separator.h"
#include "AbsBeamline/Septum.h"
#include "AbsBeamline/Solenoid.h"
#include "AbsBeamline/CyclotronValley.h"
#include "AbsBeamline/Stripper.h"
59
#include "AbsBeamline/VariableRFCavity.h"
60
#include "AbsBeamline/VariableRFCavityFringeField.h"
61

62 63
#include "AbstractObjects/Element.h"

64
#include "Beamlines/FlaggedBeamline.h"
65
#include "Elements/OpalBeamline.h"
kraus's avatar
kraus committed
66
#include "AbsBeamline/Ring.h"
gsell's avatar
gsell committed
67 68 69

#include "BeamlineGeometry/Euclid3D.h"
#include "BeamlineGeometry/PlanarArcGeometry.h"
Jianjun Yang's avatar
Jianjun Yang committed
70
#include "BeamlineGeometry/RBendGeometry.h"
ext-rogers_c's avatar
ext-rogers_c committed
71
#include "BeamlineGeometry/StraightGeometry.h"
gsell's avatar
gsell committed
72 73 74 75 76 77 78 79 80 81
#include "Beamlines/Beamline.h"

#include "Fields/BMultipoleField.h"
#include "FixedAlgebra/FTps.h"
#include "FixedAlgebra/FTpsMath.h"
#include "FixedAlgebra/FVps.h"

#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/H5PartWrapperForPC.h"
88
#include "Structure/BoundaryGeometry.h"
gsell's avatar
gsell committed
89 90 91 92 93 94 95

#include "Ctunes.h"
#include <cassert>

#include <hdf5.h>
#include "H5hut.h"

frey_m's avatar
frey_m committed
96 97
//FIXME Remove headers and dynamic_cast in readOneBunchFromFile
#include "Algorithms/PartBunch.h"
98
#ifdef ENABLE_AMR
frey_m's avatar
frey_m committed
99 100
    #include "Algorithms/AmrPartBunch.h"
#endif
frey_m's avatar
frey_m committed
101

gsell's avatar
gsell committed
102 103
class Beamline;
class PartData;
104

gsell's avatar
gsell committed
105 106 107
using Physics::pi;
using Physics::q_e;

108

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

111 112 113
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
114

snuverink_j's avatar
snuverink_j committed
115
//#define PSdim 6
gsell's avatar
gsell committed
116 117 118 119 120 121 122 123 124 125 126 127 128 129

extern Inform *gmsg;

// typedef FVector<double, PSdim> Vector;

/**
 * Constructor ParallelCyclotronTracker
 *
 * @param beamline
 * @param reference
 * @param revBeam
 * @param revTrack
 */
ParallelCyclotronTracker::ParallelCyclotronTracker(const Beamline &beamline,
130
                                                   const PartData &reference,
frey_m's avatar
frey_m committed
131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146
                                                   bool revBeam, bool revTrack)
    : Tracker(beamline, reference, revBeam, revTrack)
    , itsDataSink(nullptr)
    , itsMBDump_m(new MultiBunchDump())
    , bgf_m(nullptr)
    , lastDumpedStep_m(0)
    , eta_m(0.01)
    , myNode_m(Ippl::myNode())
    , initialLocalNum_m(0)
    , initialTotalNum_m(0)
    , onebunch_m(OpalData::getInstance()->getInputBasename() + "-onebunch.h5")
    , opalRing_m(nullptr)
    , itsStepper_mp(nullptr)
    , mode_m(MODE::UNDEFINED)
    , stepper_m(stepper::INTEGRATOR::UNDEFINED)
{
gsell's avatar
gsell committed
147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162
    itsBeamline = dynamic_cast<Beamline *>(beamline.clone());
}

/**
 * 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
163
                                                   PartBunchBase<double, 3> *bunch,
gsell's avatar
gsell committed
164 165 166
                                                   DataSink &ds,
                                                   const PartData &reference,
                                                   bool revBeam, bool revTrack,
frey_m's avatar
frey_m committed
167 168 169 170 171 172 173 174 175 176 177 178 179 180
                                                   int maxSTEPS, int timeIntegrator)
    : Tracker(beamline, bunch, reference, revBeam, revTrack)
    , itsMBDump_m(new MultiBunchDump())
    , bgf_m(nullptr)
    , maxSteps_m(maxSTEPS)
    , lastDumpedStep_m(0)
    , eta_m(0.01)
    , myNode_m(Ippl::myNode())
    , initialLocalNum_m(bunch->getLocalNum())
    , initialTotalNum_m(bunch->getTotalNum())
    , onebunch_m(OpalData::getInstance()->getInputBasename() + "-onebunch.h5")
    , opalRing_m(nullptr)
    , itsStepper_mp(nullptr)
{
gsell's avatar
gsell committed
181 182
    itsBeamline = dynamic_cast<Beamline *>(beamline.clone());
    itsDataSink = &ds;
frey_m's avatar
frey_m committed
183
    //  scaleFactor_m = itsBunch_m->getdT() * c;
gsell's avatar
gsell committed
184
    scaleFactor_m = 1;
185
    multiBunchMode_m = MB_MODE::NONE;
gsell's avatar
gsell committed
186 187 188 189 190

    IntegrationTimer_m = IpplTimings::getTimer("Integration");
    TransformTimer_m   = IpplTimings::getTimer("Frametransform");
    DumpTimer_m        = IpplTimings::getTimer("Dump");
    BinRepartTimer_m   = IpplTimings::getTimer("Binaryrepart");
191

192 193 194 195 196 197 198 199 200
    // 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;
201

202 203 204 205 206 207 208 209
    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
210 211 212 213 214 215 216
}

/**
 * Destructor ParallelCyclotronTracker
 *
 */
ParallelCyclotronTracker::~ParallelCyclotronTracker() {
217 218
    for(Component* component : myElements) {
        delete(component);
gsell's avatar
gsell committed
219
    }
220 221
    for(auto fd : FieldDimensions) {
        delete(fd);
gsell's avatar
gsell committed
222 223
    }
    delete itsBeamline;
224
    // delete opalRing_m;
gsell's avatar
gsell committed
225 226
}

227 228 229 230 231 232
/**
 * AAA
 *
 * @param none
 */
void ParallelCyclotronTracker::initializeBoundaryGeometry() {
233 234
    for(Component * component : myElements) {
        bgf_m = dynamic_cast<ElementBase *>(component)->getBoundaryGeometry();
235 236 237 238 239 240 241 242 243 244
        if(!bgf_m)
            continue;
        else
            break;
    }
    if (bgf_m) {
        itsDataSink->writeGeomToVtk(*bgf_m, std::string("data/testGeometry-00000.vtk"));
        OpalData::getInstance()->setGlobalGeometry(bgf_m);
        *gmsg << "* Boundary geometry initialized " << endl;
    }
245
}
246

247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268

/// set the working sub-mode for multi-bunch mode: "FORCE" or "AUTO"
inline
void ParallelCyclotronTracker::setMultiBunchMode(const std::string& mbmode)
{
    if ( mbmode.compare("FORCE") == 0 ) {
        *gmsg << "FORCE mode: The multi bunches will be injected consecutively "
              << "after each revolution, until get \"TURNS\" bunches." << endl;
        multiBunchMode_m = MB_MODE::FORCE;
    } else if ( mbmode.compare("AUTO") == 0 ) {
        *gmsg << "AUTO mode: The multi bunches will be injected only when "
                      << "the distance between two neighboring bunches " << endl
                      << "is below the limitation. The control parameter is set to "
                      << CoeffDBunches_m << endl;
        multiBunchMode_m = MB_MODE::AUTO;
    } else if ( mbmode.compare("NONE") == 0 )
        multiBunchMode_m = MB_MODE::NONE;
    else
        throw OpalException("ParallelCyclotronTracker::setMultiBunchMode()",
                            "MBMODE name \"" + mbmode + "\" unknown.");
}

269 270 271 272 273 274
/**
 *
 *
 * @param fn Base file name
 */
void ParallelCyclotronTracker::bgf_main_collision_test() {
275
    if(!bgf_m) return;
276

277
    Inform msg("bgf_main_collision_test ");
278

279 280 281
    /**
     *Here we check if a particles is outside the domain, flag it for deletion
     */
282

283
    Vector_t intecoords = 0.0;
284

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

288
    int triId = 0;
frey_m's avatar
frey_m committed
289 290 291
    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);
292
        if(res >= 0) {
frey_m's avatar
frey_m committed
293
            itsBunch_m->Bin[i] = -1;
294
        }
295
    }
296 297 298
}


gsell's avatar
gsell committed
299 300 301 302 303
/**
 *
 *
 * @param fn Base file name
 */
304
void ParallelCyclotronTracker::openFiles(std::string SfileName) {
gsell's avatar
gsell committed
305

306
    std::string  SfileName2 = SfileName + std::string("-Angle0.dat");
gsell's avatar
gsell committed
307
    outfTheta0_m.precision(8);
308
    outfTheta0_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
309
    outfTheta0_m.open(SfileName2.c_str());
310 311 312
    outfTheta0_m << "#  r [mm]      beta_r*gamma       "
                 << "theta [mm]      beta_theta*gamma        "
                 << "z [mm]          beta_z*gamma" << std::endl;
gsell's avatar
gsell committed
313

314
    SfileName2 = SfileName + std::string("-Angle1.dat");
gsell's avatar
gsell committed
315
    outfTheta1_m.precision(8);
316
    outfTheta1_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
317
    outfTheta1_m.open(SfileName2.c_str());
318 319 320
    outfTheta1_m << "#  r [mm]      beta_r*gamma       "
                 << "theta [mm]      beta_theta*gamma        "
                 << "z [mm]          beta_z*gamma"  << std::endl;
gsell's avatar
gsell committed
321

322
    SfileName2 = SfileName + std::string("-Angle2.dat");
gsell's avatar
gsell committed
323
    outfTheta2_m.precision(8);
324
    outfTheta2_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
325
    outfTheta2_m.open(SfileName2.c_str());
326 327 328
    outfTheta2_m << "#  r [mm]      beta_r*gamma       "
                 << "theta [mm]      beta_theta*gamma        "
                 << "z [mm]          beta_z*gamma"  << std::endl;
gsell's avatar
gsell committed
329 330 331

    // for single Particle Mode, output after each turn, to define matched initial phase ellipse.

332
    SfileName2 = SfileName + std::string("-afterEachTurn.dat");
gsell's avatar
gsell committed
333 334

    outfThetaEachTurn_m.precision(8);
335
    outfThetaEachTurn_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
336 337

    outfThetaEachTurn_m.open(SfileName2.c_str());
338 339 340
    outfTheta2_m << "#  r [mm]      beta_r*gamma       "
                 << "theta [mm]      beta_theta*gamma        "
                 << "z [mm]          beta_z*gamma"  << std::endl;
gsell's avatar
gsell committed
341 342 343 344 345 346 347 348 349 350 351 352 353 354 355
}

/**
 * Close all files related to
 * special output in the Cyclotron
 * mode.
 */
void ParallelCyclotronTracker::closeFiles() {

    outfTheta0_m.close();
    outfTheta1_m.close();
    outfTheta2_m.close();
    outfThetaEachTurn_m.close();
}

356
/**
357 358 359
 *
 * @param ring
 */
kraus's avatar
kraus committed
360
void ParallelCyclotronTracker::visitRing(const Ring &ring) {
Daniel Winklehner's avatar
Daniel Winklehner committed
361

362
    *gmsg << "* ----------------------------- Adding Ring ------------------------------ *" << endl;
Daniel Winklehner's avatar
Daniel Winklehner committed
363

snuverink_j's avatar
snuverink_j committed
364
    delete opalRing_m;
Daniel Winklehner's avatar
Daniel Winklehner committed
365

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

368
    myElements.push_back(opalRing_m);
Daniel Winklehner's avatar
Daniel Winklehner committed
369

frey_m's avatar
frey_m committed
370
    opalRing_m->initialise(itsBunch_m);
371 372 373 374

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

376
    if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
kraus's avatar
kraus committed
377
        throw OpalException("Error in ParallelCyclotronTracker::visitRing",
378 379
                            "PHIINIT is out of [-180, 180)!");
    }
Daniel Winklehner's avatar
Daniel Winklehner committed
380 381

    referenceZ = 0.0;
382
    referencePz = 0.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
383 384 385 386

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

387 388
    if(referencePtot < 0.0)
        referencePt *= -1.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
389

390 391
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
gsell's avatar
gsell committed
392

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

395
    buildupFieldList(BcParameter, ElementBase::RING, opalRing_m);
396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411

    // 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
412 413 414 415 416 417 418 419

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

420
    *gmsg << "* -------------------------- Adding Cyclotron ---------------------------- *" << endl;
gsell's avatar
gsell committed
421

422 423
    Cyclotron *elptr = dynamic_cast<Cyclotron *>(cycl.clone());
    myElements.push_back(elptr);
424

425
    // Is this a Spiral Inflector Simulation? If yes, we'll give the user some
426 427 428 429 430 431
    // useful information
    spiral_flag = elptr->getSpiralFlag();

    if(spiral_flag) {

        *gmsg << endl << "* This is a Spiral Inflector Simulation! This means the following:" << endl;
432
        *gmsg         << "* 1.) It is up to the user to provide appropriate geometry, electric and magnetic fields!" << endl;
433 434
        *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
435
        *gmsg         << "* 2.) For high currents it is strongly recommended to use the SAAMG fieldsolver," << endl;
436 437
        *gmsg         << "*     FFT does not give the correct results (boundaty conditions are missing)." << endl;
        *gmsg         << "* 3.) The whole geometry will be meshed and used for the fieldsolve." << endl;
438 439
        *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;
440 441 442 443 444 445 446
        *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;
        numBunch_m = 1;

    }

447
    // Fresh run (no restart):
448
    if(!OpalData::getInstance()->inRestartRun()) {
449

450
        // Get reference values from cyclotron element
451
        // For now, these are still stored in mm. should be the only ones. -DW
452 453 454 455
        referenceR     = elptr->getRinit();
        referenceTheta = elptr->getPHIinit();
        referenceZ     = elptr->getZinit();
        referencePr    = elptr->getPRinit();
456 457
        referencePz    = elptr->getPZinit();

458
        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
459 460
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                "PHIINIT is out of [-180, 180)!");
461 462 463 464 465
        }

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

        // Calculate reference azimuthal (tangential) momentum from total-, z- and radial momentum:
466
        float insqrt = referencePtot * referencePtot - \
467
            referencePr * referencePr - referencePz * referencePz;
468 469 470 471 472

        if(insqrt < 0) {

            if(insqrt > -1.0e-10) {

473
                referencePt = 0.0;
474 475 476

            } else {

477
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
478
                                    "Pt imaginary!");
479 480 481 482 483 484 485
            }

        } else {

            referencePt = sqrt(insqrt);
        }

486
        if(referencePtot < 0.0)
Daniel Winklehner's avatar
Daniel Winklehner committed
487
            referencePt *= -1.0;
488 489
        // End calculate referencePt

490
        // Restart a run:
491 492
    } else {

493
        // If the user wants to save the restarted run in local frame,
494
        // make sure the previous h5 file was local too
495
      if (Options::psDumpFrame != Options::GLOBAL) {
winklehner_d's avatar
-DW  
winklehner_d committed
496
        if (!previousH5Local) {
497 498
                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
499
        }
500 501
            // 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
502 503
      } else {
        if (previousH5Local) {
504 505
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                    "You are trying a global restart from a local h5 file!");
506
            }
507
        }
508

509
        // Adjust some of the reference variables from the h5 file
510 511
        referencePhi *= Physics::deg2rad;
        referencePsi *= Physics::deg2rad;
512
        referencePtot = bega;
513
        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
514
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
515 516
                                "PHIINIT is out of [-180, 180)!");
        }
517 518
    }

519 520
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
521

522
    *gmsg << endl;
adelmann's avatar
adelmann committed
523
    *gmsg << "* Bunch global starting position:" << endl;
524 525
    *gmsg << "* RINIT = " << referenceR  << " [mm]" << endl;
    *gmsg << "* PHIINIT = " << referenceTheta << " [deg]" << endl;
526
    *gmsg << "* ZINIT = " << referenceZ << " [mm]" << endl;
527
    *gmsg << endl;
adelmann's avatar
adelmann committed
528
    *gmsg << "* Bunch global starting momenta:" << endl;
gsell's avatar
gsell committed
529 530
    *gmsg << "* Initial gamma = " << itsReference.getGamma() << endl;
    *gmsg << "* Initial beta = " << itsReference.getBeta() << endl;
531
    *gmsg << "* Reference total momentum (beta * gamma) = " << referencePtot * 1000.0 << " [MCU]" << endl;
532 533 534
    *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;
535
    *gmsg << endl;
adelmann's avatar
adelmann committed
536

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

540 541 542
    // 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
543

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

547
    std::string type = elptr->getCyclotronType();
548
    *gmsg << "* Type of cyclotron = " << type << " " << endl;
549

550 551
    double rmin = elptr->getMinR();
    double rmax = elptr->getMaxR();
adelmann's avatar
adelmann committed
552
    *gmsg << "* Radial aperture = " << rmin << " ... " << rmax<<" [m] "<< endl;
553 554 555

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

    double h = elptr->getCyclHarm();
559
    *gmsg << "* Number of trimcoils = " << elptr->getNumberOfTrimcoils() << endl;
560
    *gmsg << "* Harmonic number h = " << h << " " << endl;
561

gsell's avatar
gsell committed
562 563 564 565 566 567 568 569 570
    /**
     * 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
     * fieldflag = 5, readin FFAG format file for MSU/FNAL FFAG
     * fieldflag = 6, readin both median plane B field map and 3D E field map of RF cavity for compact cyclotron
571
     * fieldflag = 7, read in fields for Daniel's synchrocyclotron simulations
gsell's avatar
gsell committed
572
     */
573
    int fieldflag = elptr->getFieldFlag(type);
gsell's avatar
gsell committed
574

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

snuverink_j's avatar
snuverink_j committed
578
    double BcParameter[8] = {};
579

580 581
    BcParameter[0] = 0.001 * elptr->getRmin();
    BcParameter[1] = 0.001 * elptr->getRmax();
gsell's avatar
gsell committed
582

583
    // Store inner radius and outer radius of cyclotron field map in the list
584
    buildupFieldList(BcParameter, ElementBase::CYCLOTRON, elptr);
gsell's avatar
gsell committed
585 586 587 588 589 590 591 592 593 594 595 596 597 598 599
}

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

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

602
    *gmsg << "* --------- Collimator -----------------------------" << endl;
gsell's avatar
gsell committed
603

604
    CCollimator* elptr = dynamic_cast<CCollimator *>(coll.clone());
605
    myElements.push_back(elptr);
gsell's avatar
gsell committed
606

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

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

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

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

619
    double zstart = elptr->getZStart();
adelmann's avatar
adelmann committed
620
    *gmsg << "* Zstart= " << zstart << " [mm]" << endl;
621 622

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

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

frey_m's avatar
frey_m committed
628
    elptr->initialise(itsBunch_m);
gsell's avatar
gsell committed
629

snuverink_j's avatar
snuverink_j committed
630
    double BcParameter[8] = {};
631

632 633 634 635 636 637
    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 ;

638
    buildupFieldList(BcParameter, ElementBase::CCOLLIMATOR, elptr);
gsell's avatar
gsell committed
639 640 641 642 643 644 645 646 647 648 649 650
}

/**
 *
 *
 * @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
651 652 653 654 655 656 657 658 659 660 661 662
/**
 *
 *
 * @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
663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682
/**
 *
 *
 * @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()));
}

683 684
/**
 *
685 686
 *
 *  @param
687 688 689 690 691
 */
void ParallelCyclotronTracker::visitFlexibleCollimator(const FlexibleCollimator &) {

}

gsell's avatar
gsell committed
692 693 694 695 696 697 698 699 700 701
/**
 *
 *
 * @param lamb
 */
void ParallelCyclotronTracker::visitLambertson(const Lambertson &lamb) {
    *gmsg << "In Lambertson; L= " << lamb.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Lambertson *>(lamb.clone()));
}

702 703 704
void ParallelCyclotronTracker::visitOffset(const Offset & off) {
    if (opalRing_m == NULL)
        throw OpalException(
705 706
                            "ParallelCylcotronTracker::visitOffset",
                            "Attempt to place an offset when Ring not defined");
707 708
    Offset* offNonConst = const_cast<Offset*>(&off);
    offNonConst->updateGeometry(opalRing_m->getNextPosition(),
709
                                opalRing_m->getNextNormal());
710 711 712
    opalRing_m->appendElement(off);
}

gsell's avatar
gsell committed
713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744
/**
 *
 *
 * @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
745 746 747 748 749 750 751
/**
 *
 *
 * @param multT
 */
void ParallelCyclotronTracker::visitMultipoleT(const MultipoleT &multT) {
    *gmsg << "Adding MultipoleT" << endl;
752
    if (opalRing_m != NULL) {
ext-rogers_c's avatar
ext-rogers_c committed
753
        opalRing_m->appendElement(multT);
754
    } else {
ext-rogers_c's avatar
ext-rogers_c committed
755 756
        throw OpalException("ParallelCyclotronTracker::visitMultipoleT",
                            "Need to define a RINGDEFINITION to use MultipoleT element");
757
    }
ext-rogers_c's avatar
ext-rogers_c committed
758 759 760
    myElements.push_back(dynamic_cast<MultipoleT *>(multT.clone()));
}

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 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808
/**
 *
 *
 * @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
809 810 811 812 813 814
/**
 *
 *
 * @param prob
 */
void ParallelCyclotronTracker::visitProbe(const Probe &prob) {
815
    *gmsg << "* -----------  Probe -------------------------------" << endl;
816 817
    Probe *elptr = dynamic_cast<Probe *>(prob.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
818

819
    double xstart = elptr->getXstart();
820
    *gmsg << "XStart= " << xstart << " [mm]" << endl;
gsell's avatar
gsell committed
821

822
    double xend = elptr->getXend();
823
    *gmsg << "XEnd= " << xend << " [mm]" << endl;
gsell's avatar
gsell committed
824

825
    double ystart = elptr->getYstart();
826
    *gmsg << "YStart= " << ystart << " [mm]" << endl;
gsell's avatar
gsell committed
827

828
    double yend = elptr->getYend();
829
    *gmsg << "YEnd= " << yend << " [mm]" << endl;
gsell's avatar
gsell committed
830 831

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

snuverink_j's avatar
snuverink_j committed
834
    double BcParameter[8] = {};
835

836 837 838 839
    BcParameter[0] = 0.001 * xstart ;
    BcParameter[1] = 0.001 * xend;
    BcParameter[2] = 0.001 * ystart ;
    BcParameter[3] = 0.001 * yend;
840
    BcParameter[4] = 0.001 ; // width
gsell's avatar
gsell committed
841 842

    // store probe parameters in the list
843
    buildupFieldList(BcParameter, ElementBase::PROBE, elptr);
gsell's avatar
gsell committed
844 845 846 847 848 849 850 851 852 853 854 855
}

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

856 857 858 859 860 861
void ParallelCyclotronTracker::visitSBend3D(const SBend3D &bend) {
    *gmsg << "Adding SBend3D" << endl;
    if (opalRing_m != NULL)
        opalRing_m->appendElement(bend);
    else
        throw OpalException("ParallelCyclotronTracker::visitSBend3D",
862
                            "Need to define a RINGDEFINITION to use SBend3D element");
863 864
}

865 866
void ParallelCyclotronTracker::visitScalingFFAGMagnet(const ScalingFFAGMagnet &bend) {
    *gmsg << "Adding ScalingFFAGMagnet" << endl;
ext-rogers_c's avatar
ext-rogers_c committed
867
    if (opalRing_m != NULL) {
868
        opalRing_m->appendElement(bend);
ext-rogers_c's avatar
ext-rogers_c committed
869
    } else {
870 871
        throw OpalException("ParallelCyclotronTracker::visitScalingFFAGMagnet",
                            "Need to define a RINGDEFINITION to use ScalingFFAGMagnet element");
ext-rogers_c's avatar
ext-rogers_c committed
872
    }
873 874
}

875 876 877 878 879 880
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",
881
                            "Need to define a RINGDEFINITION to use VariableRFCavity element");
882 883 884 885 886 887 888 889 890 891
}

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");
892 893
}

gsell's avatar
gsell committed
894 895 896 897 898 899 900 901
/**
 *
 *
 * @param as
 */
void ParallelCyclotronTracker::visitRFCavity(const RFCavity &as) {

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

903 904
    RFCavity *elptr = dynamic_cast<RFCavity *>(as.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
905 906 907 908 909 910 911

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

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

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

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

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

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

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

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

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

936 937 938
    /*
      Setup time dependence and in case of no
      timedependence use a polynom with  a_0 = 1 and a_k = 0, k = 1,2,3.
939
    */
940

941 942 943 944
    std::shared_ptr<AbstractTimeDependence> freq_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> ampl_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> phase_atd = nullptr;

945
    dvector_t  unityVec;
946 947 948 949
    unityVec.push_back(1.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
950

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

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