ParallelCyclotronTracker.cpp 120 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>
25 26

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

29
#include "AbsBeamline/CCollimator.h"
gsell's avatar
gsell committed
30 31
#include "AbsBeamline/Corrector.h"
#include "AbsBeamline/Cyclotron.h"
32
#include "AbsBeamline/CyclotronValley.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"
47
#include "AbsBeamline/PluginElement.h"
gsell's avatar
gsell committed
48 49
#include "AbsBeamline/Probe.h"
#include "AbsBeamline/RBend.h"
50
#include "AbsBeamline/Ring.h"
gsell's avatar
gsell committed
51 52 53
#include "AbsBeamline/RFCavity.h"
#include "AbsBeamline/RFQuadrupole.h"
#include "AbsBeamline/SBend.h"
54
#include "AbsBeamline/SBend3D.h"
ext-rogers_c's avatar
ext-rogers_c committed
55
#include "AbsBeamline/ScalingFFAMagnet.h"
gsell's avatar
gsell committed
56 57 58 59
#include "AbsBeamline/Separator.h"
#include "AbsBeamline/Septum.h"
#include "AbsBeamline/Solenoid.h"
#include "AbsBeamline/Stripper.h"
60
#include "AbsBeamline/VariableRFCavity.h"
61
#include "AbsBeamline/VariableRFCavityFringeField.h"
62

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

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

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

#include "Physics/Physics.h"

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

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

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

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

85

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

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

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

130 131 132
    if ( numBunch > 1 ) {
        mbHandler_m = std::unique_ptr<MultiBunchHandler>(
            new MultiBunchHandler(bunch, numBunch, mbEta,
133
                                  mbPara, mbMode, mbBinning, ds)
134 135 136
        );
    }

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

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

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

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

frey_m's avatar
frey_m committed
179

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

188
    Inform msg("bgf_main_collision_test ");
189

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

194
    Vector_t intecoords = 0.0;
195

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

199
    int triId = 0;
frey_m's avatar
frey_m committed
200 201 202
    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);
203
        if(res >= 0) {
frey_m's avatar
frey_m committed
204
            itsBunch_m->Bin[i] = -1;
205
        }
206
    }
207 208
}

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


230 231 232 233
double ParallelCyclotronTracker::computePathLengthUpdate(const double& dt,
                                                         short bunchNr)
{
    double dotP = 0.0;
234
    if ( Options::psDumpFrame == Options::BUNCH_MEAN || isMultiBunch()) {
235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258
        for(unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) {
            // take all particles if bunchNr <= -1
            if ( bunchNr > -1 && itsBunch_m->bunchNum[i] != bunchNr)
                continue;

            dotP += dot(itsBunch_m->P[i], itsBunch_m->P[i]);
        }

        allreduce(dotP, 1, std::plus<double>());

        dotP /= itsBunch_m->getTotalNum();

    } else if ( itsBunch_m->getLocalNum() == 0 ) {
        return 0.0;
    } else {
        dotP = dot(itsBunch_m->P[0], itsBunch_m->P[0]);
    }

    double const gamma = std::sqrt(1.0 + dotP);
    double const beta  = std::sqrt(dotP) / gamma;
    return c_mmtns * dt * 1.0e-3 * beta; // unit: m
}


259

gsell's avatar
gsell committed
260 261 262 263 264
/**
 *
 *
 * @param fn Base file name
 */
265
void ParallelCyclotronTracker::openFiles(size_t numFiles, std::string SfileName) {
gsell's avatar
gsell committed
266

267
    for (unsigned int i=0; i<numFiles; i++) {
268 269 270 271 272 273 274 275
        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
276

277 278 279 280 281 282
        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;
283
    }
gsell's avatar
gsell committed
284 285 286 287 288 289 290 291
}

/**
 * Close all files related to
 * special output in the Cyclotron
 * mode.
 */
void ParallelCyclotronTracker::closeFiles() {
292
    for (auto & file : outfTheta_m) {
293
        file->close();
294
    }
gsell's avatar
gsell committed
295 296
}

297
/**
298 299 300
 *
 * @param ring
 */
kraus's avatar
kraus committed
301
void ParallelCyclotronTracker::visitRing(const Ring &ring) {
Daniel Winklehner's avatar
Daniel Winklehner committed
302

303
    *gmsg << "* ----------------------------- Adding Ring ------------------------------ *" << endl;
Daniel Winklehner's avatar
Daniel Winklehner committed
304

snuverink_j's avatar
snuverink_j committed
305
    delete opalRing_m;
Daniel Winklehner's avatar
Daniel Winklehner committed
306

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

309
    myElements.push_back(opalRing_m);
Daniel Winklehner's avatar
Daniel Winklehner committed
310

frey_m's avatar
frey_m committed
311
    opalRing_m->initialise(itsBunch_m);
312 313 314 315

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

317
    if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
kraus's avatar
kraus committed
318
        throw OpalException("Error in ParallelCyclotronTracker::visitRing",
319 320
                            "PHIINIT is out of [-180, 180)!");
    }
Daniel Winklehner's avatar
Daniel Winklehner committed
321 322

    referenceZ = 0.0;
323
    referencePz = 0.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
324 325 326 327

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

328 329
    if(referencePtot < 0.0)
        referencePt *= -1.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
330

331 332
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
gsell's avatar
gsell committed
333

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

336
    buildupFieldList(BcParameter, ElementBase::RING, opalRing_m);
337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352

    // 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
353 354 355 356 357 358 359 360

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

361
    *gmsg << "* -------------------------- Adding Cyclotron ---------------------------- *" << endl;
gsell's avatar
gsell committed
362

363 364
    Cyclotron *elptr = dynamic_cast<Cyclotron *>(cycl.clone());
    myElements.push_back(elptr);
365

366
    // Is this a Spiral Inflector Simulation? If yes, we'll give the user some
367 368 369 370 371 372
    // useful information
    spiral_flag = elptr->getSpiralFlag();

    if(spiral_flag) {

        *gmsg << endl << "* This is a Spiral Inflector Simulation! This means the following:" << endl;
373
        *gmsg         << "* 1.) It is up to the user to provide appropriate geometry, electric and magnetic fields!" << endl;
374 375
        *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
376
        *gmsg         << "* 2.) For high currents it is strongly recommended to use the SAAMG fieldsolver," << endl;
377 378
        *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;
379 380
        *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;
381 382 383
        *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;
384
        mbHandler_m->setTotalNumBunch(1);
385 386
    }

387
    // Fresh run (no restart):
388
    if(!OpalData::getInstance()->inRestartRun()) {
389

390
        // Get reference values from cyclotron element
391
        // For now, these are still stored in mm. should be the only ones. -DW
392 393 394 395
        referenceR     = elptr->getRinit();
        referenceTheta = elptr->getPHIinit();
        referenceZ     = elptr->getZinit();
        referencePr    = elptr->getPRinit();
396 397
        referencePz    = elptr->getPZinit();

398
        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
399 400
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                "PHIINIT is out of [-180, 180)!");
401 402 403 404 405
        }

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

        // Calculate reference azimuthal (tangential) momentum from total-, z- and radial momentum:
406
        float insqrt = referencePtot * referencePtot - \
407
            referencePr * referencePr - referencePz * referencePz;
408 409 410 411 412

        if(insqrt < 0) {

            if(insqrt > -1.0e-10) {

413
                referencePt = 0.0;
414 415 416

            } else {

417
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
418
                                    "Pt imaginary!");
419 420 421 422 423 424 425
            }

        } else {

            referencePt = sqrt(insqrt);
        }

426
        if(referencePtot < 0.0)
Daniel Winklehner's avatar
Daniel Winklehner committed
427
            referencePt *= -1.0;
428 429
        // End calculate referencePt

430
        // Restart a run:
431 432
    } else {

433
        // If the user wants to save the restarted run in local frame,
434
        // make sure the previous h5 file was local too
435
      if (Options::psDumpFrame != Options::GLOBAL) {
winklehner_d's avatar
-DW  
winklehner_d committed
436
        if (!previousH5Local) {
437 438
                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
439
        }
440 441
            // 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
442 443
      } else {
        if (previousH5Local) {
444 445
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                    "You are trying a global restart from a local h5 file!");
446
            }
447
        }
448

449
        // Adjust some of the reference variables from the h5 file
450 451
        referencePhi *= Physics::deg2rad;
        referencePsi *= Physics::deg2rad;
452
        referencePtot = bega;
453
        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
454
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
455 456
                                "PHIINIT is out of [-180, 180)!");
        }
457 458
    }

459 460
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
461

462
    *gmsg << endl;
adelmann's avatar
adelmann committed
463
    *gmsg << "* Bunch global starting position:" << endl;
464 465
    *gmsg << "* RINIT = " << referenceR  << " [mm]" << endl;
    *gmsg << "* PHIINIT = " << referenceTheta << " [deg]" << endl;
466
    *gmsg << "* ZINIT = " << referenceZ << " [mm]" << endl;
467
    *gmsg << endl;
adelmann's avatar
adelmann committed
468
    *gmsg << "* Bunch global starting momenta:" << endl;
gsell's avatar
gsell committed
469 470
    *gmsg << "* Initial gamma = " << itsReference.getGamma() << endl;
    *gmsg << "* Initial beta = " << itsReference.getBeta() << endl;
471
    *gmsg << "* Reference total momentum (beta * gamma) = " << referencePtot * 1000.0 << " [MCU]" << endl;
472 473 474
    *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;
475
    *gmsg << endl;
adelmann's avatar
adelmann committed
476

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

480 481 482
    // 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
483

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

487
    std::string type = elptr->getCyclotronType();
488
    *gmsg << "* Type of cyclotron = " << type << " " << endl;
489

490 491
    double rmin = elptr->getMinR();
    double rmax = elptr->getMaxR();
adelmann's avatar
adelmann committed
492
    *gmsg << "* Radial aperture = " << rmin << " ... " << rmax<<" [m] "<< endl;
493 494 495

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

    double h = elptr->getCyclHarm();
499
    *gmsg << "* Number of trimcoils = " << elptr->getNumberOfTrimcoils() << endl;
500
    *gmsg << "* Harmonic number h = " << h << " " << endl;
501

gsell's avatar
gsell committed
502 503 504 505 506 507 508
    /**
     * 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
509
     * fieldflag = 5, readin FFA format file for MSU/FNAL FFA
gsell's avatar
gsell committed
510
     * fieldflag = 6, readin both median plane B field map and 3D E field map of RF cavity for compact cyclotron
511
     * fieldflag = 7, read in fields for Daniel's synchrocyclotron simulations
gsell's avatar
gsell committed
512
     */
513
    int fieldflag = elptr->getFieldFlag(type);
gsell's avatar
gsell committed
514

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

snuverink_j's avatar
snuverink_j committed
518
    double BcParameter[8] = {};
519

520 521
    BcParameter[0] = 0.001 * elptr->getRmin();
    BcParameter[1] = 0.001 * elptr->getRmax();
gsell's avatar
gsell committed
522

523
    // Store inner radius and outer radius of cyclotron field map in the list
524
    buildupFieldList(BcParameter, ElementBase::CYCLOTRON, elptr);
gsell's avatar
gsell committed
525 526 527 528 529 530 531 532 533 534 535 536 537 538 539
}

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

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

542
    *gmsg << "* --------- Collimator -----------------------------" << endl;
gsell's avatar
gsell committed
543

544
    CCollimator* elptr = dynamic_cast<CCollimator *>(coll.clone());
545
    myElements.push_back(elptr);
gsell's avatar
gsell committed
546

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

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

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

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

559
    double zstart = elptr->getZStart();
adelmann's avatar
adelmann committed
560
    *gmsg << "* Zstart= " << zstart << " [mm]" << endl;
561 562

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

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

frey_m's avatar
frey_m committed
568
    elptr->initialise(itsBunch_m);
gsell's avatar
gsell committed
569

snuverink_j's avatar
snuverink_j committed
570
    double BcParameter[8] = {};
571

572 573 574 575 576 577
    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 ;

578
    buildupFieldList(BcParameter, ElementBase::CCOLLIMATOR, elptr);
gsell's avatar
gsell committed
579 580 581 582 583 584 585 586 587 588 589 590
}

/**
 *
 *
 * @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
591 592 593 594 595 596 597 598 599 600 601 602
/**
 *
 *
 * @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
603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622
/**
 *
 *
 * @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()));
}

623 624
/**
 *
625 626
 *
 *  @param
627 628 629 630 631
 */
void ParallelCyclotronTracker::visitFlexibleCollimator(const FlexibleCollimator &) {

}

gsell's avatar
gsell committed
632 633 634 635 636 637 638 639 640 641
/**
 *
 *
 * @param lamb
 */
void ParallelCyclotronTracker::visitLambertson(const Lambertson &lamb) {
    *gmsg << "In Lambertson; L= " << lamb.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Lambertson *>(lamb.clone()));
}

642 643 644
void ParallelCyclotronTracker::visitOffset(const Offset & off) {
    if (opalRing_m == NULL)
        throw OpalException(
645 646
                            "ParallelCylcotronTracker::visitOffset",
                            "Attempt to place an offset when Ring not defined");
647 648
    Offset* offNonConst = const_cast<Offset*>(&off);
    offNonConst->updateGeometry(opalRing_m->getNextPosition(),
649
                                opalRing_m->getNextNormal());
650 651 652
    opalRing_m->appendElement(off);
}

gsell's avatar
gsell committed
653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684
/**
 *
 *
 * @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
685 686 687 688 689 690 691
/**
 *
 *
 * @param multT
 */
void ParallelCyclotronTracker::visitMultipoleT(const MultipoleT &multT) {
    *gmsg << "Adding MultipoleT" << endl;
692
    if (opalRing_m != NULL) {
ext-rogers_c's avatar
ext-rogers_c committed
693
        opalRing_m->appendElement(multT);
694
    } else {
ext-rogers_c's avatar
ext-rogers_c committed
695 696
        throw OpalException("ParallelCyclotronTracker::visitMultipoleT",
                            "Need to define a RINGDEFINITION to use MultipoleT element");
697
    }
ext-rogers_c's avatar
ext-rogers_c committed
698 699 700
    myElements.push_back(dynamic_cast<MultipoleT *>(multT.clone()));
}

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 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748
/**
 *
 *
 * @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
749 750 751 752 753 754
/**
 *
 *
 * @param prob
 */
void ParallelCyclotronTracker::visitProbe(const Probe &prob) {
755
    *gmsg << "* -----------  Probe -------------------------------" << endl;
756 757
    Probe *elptr = dynamic_cast<Probe *>(prob.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
758

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

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

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

768
    double yend = elptr->getYEnd();
769
    *gmsg << "YEnd= " << yend << " [mm]" << endl;
gsell's avatar
gsell committed
770 771

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

snuverink_j's avatar
snuverink_j committed
774
    double BcParameter[8] = {};
775

776 777 778 779
    BcParameter[0] = 0.001 * xstart ;
    BcParameter[1] = 0.001 * xend;
    BcParameter[2] = 0.001 * ystart ;
    BcParameter[3] = 0.001 * yend;
780
    BcParameter[4] = 0.001 ; // width
gsell's avatar
gsell committed
781 782

    // store probe parameters in the list
783
    buildupFieldList(BcParameter, ElementBase::PROBE, elptr);
gsell's avatar
gsell committed
784 785 786 787 788 789 790 791 792 793 794 795
}

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

796 797 798 799 800 801
void ParallelCyclotronTracker::visitSBend3D(const SBend3D &bend) {
    *gmsg << "Adding SBend3D" << endl;
    if (opalRing_m != NULL)
        opalRing_m->appendElement(bend);
    else
        throw OpalException("ParallelCyclotronTracker::visitSBend3D",
802
                            "Need to define a RINGDEFINITION to use SBend3D element");
803 804
}

ext-rogers_c's avatar
ext-rogers_c committed
805 806
void ParallelCyclotronTracker::visitScalingFFAMagnet(const ScalingFFAMagnet &bend) {
    *gmsg << "Adding ScalingFFAMagnet" << endl;
ext-rogers_c's avatar
ext-rogers_c committed
807
    if (opalRing_m != NULL) {
808
        opalRing_m->appendElement(bend);
ext-rogers_c's avatar
ext-rogers_c committed
809
    } else {
ext-rogers_c's avatar
ext-rogers_c committed
810 811
        throw OpalException("ParallelCyclotronTracker::visitScalingFFAMagnet",
                            "Need to define a RINGDEFINITION to use ScalingFFAMagnet element");
ext-rogers_c's avatar
ext-rogers_c committed
812
    }
813 814
}

815 816 817 818 819 820
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",
821
                            "Need to define a RINGDEFINITION to use VariableRFCavity element");
822 823 824 825 826 827 828 829 830 831
}

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

gsell's avatar
gsell committed
834 835 836 837 838 839 840 841
/**
 *
 *
 * @param as
 */
void ParallelCyclotronTracker::visitRFCavity(const RFCavity &as) {

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

843 844
    RFCavity *elptr = dynamic_cast<RFCavity *>(as.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
845 846 847 848 849 850 851

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

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

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

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

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

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

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

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

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

876 877 878
    /*
      Setup time dependence and in case of no
      timedependence use a polynom with  a_0 = 1 and a_k = 0, k = 1,2,3.
879
    */
880

881 882 883 884
    std::shared_ptr<AbstractTimeDependence> freq_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> ampl_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> phase_atd = nullptr;

885
    dvector_t  unityVec;
886 887 888 889
    unityVec.push_back(1.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
890

891
    if (elptr->getFrequencyModelName() != "") {
892 893
        freq_atd = AbstractTimeDependence::getTimeDependence(elptr->getFrequencyModelName());
        *gmsg << "* Variable frequency RF Model name " << elptr->getFrequencyModelName() << endl;
894
    }
895 896
    else
        freq_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
897 898

    if (elptr->getAmplitudeModelName() != "") {
899 900
        ampl_atd = AbstractTimeDependence::getTimeDependence(elptr->getAmplitudeModelName());
        *gmsg << "* Variable amplitude RF Model name " << elptr->getAmplitudeModelName() << endl;
901
    }
902 903
    else
        ampl_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
904 905

    if (elptr->getPhaseModelName() != "") {
906 907
        phase_atd = AbstractTimeDependence::getTimeDependence(elptr->getPhaseModelName());
        *gmsg << "* Variable phase RF Model name " << elptr->getPhaseModelName() << endl;
908
    }
909 910
    else
        phase_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
911

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

snuverink_j's avatar
snuverink_j committed
915
    double BcParameter[8] = {};
916

917 918
    BcParameter[0] = 0.001 * rmin;
    BcParameter[1] = 0.001 * rmax;
919
    BcParameter[2] = 0.001 * pdis;
gsell's avatar
gsell committed
920