ParallelCyclotronTracker.cpp 118 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
    : Tracker(beamline, bunch, reference, revBeam, revTrack)
    , bgf_m(nullptr)
    , maxSteps_m(maxSTEPS)
120 121
    , mbHandler_m(new MultiBunchHandler(bunch, numBunch, mbEta,
                                        mbPara, mbMode, mbBinning))
frey_m's avatar
frey_m committed
122 123 124 125 126 127 128
    , 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
129 130 131 132 133 134 135
    itsBeamline = dynamic_cast<Beamline *>(beamline.clone());
    itsDataSink = &ds;

    IntegrationTimer_m = IpplTimings::getTimer("Integration");
    TransformTimer_m   = IpplTimings::getTimer("Frametransform");
    DumpTimer_m        = IpplTimings::getTimer("Dump");
    BinRepartTimer_m   = IpplTimings::getTimer("Binaryrepart");
136 137
    PluginElemTimer_m  = IpplTimings::getTimer("PluginElements");
    DelParticleTimer_m = IpplTimings::getTimer("DeleteParticles");
138

139 140 141 142 143 144 145 146 147
    // 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;
148

149 150 151 152 153 154 155 156
    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
157 158 159 160 161 162 163
}

/**
 * Destructor ParallelCyclotronTracker
 *
 */
ParallelCyclotronTracker::~ParallelCyclotronTracker() {
164 165
    for(Component* component : myElements) {
        delete(component);
gsell's avatar
gsell committed
166
    }
167 168
    for(auto fd : FieldDimensions) {
        delete(fd);
gsell's avatar
gsell committed
169 170
    }
    delete itsBeamline;
171
    // delete opalRing_m;
gsell's avatar
gsell committed
172 173
}

frey_m's avatar
frey_m committed
174

175 176 177 178 179 180
/**
 *
 *
 * @param fn Base file name
 */
void ParallelCyclotronTracker::bgf_main_collision_test() {
181
    if(!bgf_m) return;
182

183
    Inform msg("bgf_main_collision_test ");
184

185
    /**
186
     *Here we check if a particle is outside the domain, flag it for deletion
187
     */
188

189
    Vector_t intecoords = 0.0;
190

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

194
    int triId = 0;
frey_m's avatar
frey_m committed
195 196 197
    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);
198
        if(res >= 0) {
frey_m's avatar
frey_m committed
199
            itsBunch_m->Bin[i] = -1;
200
        }
201
    }
202 203
}

204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221
// only used for dumping into stat file
void ParallelCyclotronTracker::dumpAngle_m(const double& theta) {
    if ( prevAzimuth_m < 0.0 ) { // only at first occurrence
        azimuth_m = theta;
    } else {
        double dtheta = theta - prevAzimuth_m;
        if ( dtheta < 0 ) {
            dtheta += 360.0;
        }
        if ( dtheta > 180 ) { // rotating clockwise, reduce angle
            dtheta -= 360;
        }
        azimuth_m += dtheta;
    }
    prevAzimuth_m = theta;
}


222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250
double ParallelCyclotronTracker::computePathLengthUpdate(const double& dt,
                                                         short bunchNr)
{
    double dotP = 0.0;
    if ( Options::psDumpFrame == Options::BUNCH_MEAN || mbHandler_m->isMultiBunch()) {
        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
}


251

gsell's avatar
gsell committed
252 253 254 255 256
/**
 *
 *
 * @param fn Base file name
 */
257
void ParallelCyclotronTracker::openFiles(size_t numFiles, std::string SfileName) {
gsell's avatar
gsell committed
258

259
    for (unsigned int i=0; i<numFiles; i++) {
260 261 262 263 264 265 266 267
        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
268

269 270 271 272 273 274
        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;
275
    }
gsell's avatar
gsell committed
276 277 278 279 280 281 282 283
}

/**
 * Close all files related to
 * special output in the Cyclotron
 * mode.
 */
void ParallelCyclotronTracker::closeFiles() {
284
    for (auto & file : outfTheta_m) {
285
        file->close();
286
    }
gsell's avatar
gsell committed
287 288
}

289
/**
290 291 292
 *
 * @param ring
 */
kraus's avatar
kraus committed
293
void ParallelCyclotronTracker::visitRing(const Ring &ring) {
Daniel Winklehner's avatar
Daniel Winklehner committed
294

295
    *gmsg << "* ----------------------------- Adding Ring ------------------------------ *" << endl;
Daniel Winklehner's avatar
Daniel Winklehner committed
296

snuverink_j's avatar
snuverink_j committed
297
    delete opalRing_m;
Daniel Winklehner's avatar
Daniel Winklehner committed
298

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

301
    myElements.push_back(opalRing_m);
Daniel Winklehner's avatar
Daniel Winklehner committed
302

frey_m's avatar
frey_m committed
303
    opalRing_m->initialise(itsBunch_m);
304 305 306 307

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

309
    if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
kraus's avatar
kraus committed
310
        throw OpalException("Error in ParallelCyclotronTracker::visitRing",
311 312
                            "PHIINIT is out of [-180, 180)!");
    }
Daniel Winklehner's avatar
Daniel Winklehner committed
313 314

    referenceZ = 0.0;
315
    referencePz = 0.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
316 317 318 319

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

320 321
    if(referencePtot < 0.0)
        referencePt *= -1.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
322

323 324
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
gsell's avatar
gsell committed
325

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

328
    buildupFieldList(BcParameter, ElementBase::RING, opalRing_m);
329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344

    // 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
345 346 347 348 349 350 351 352

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

353
    *gmsg << "* -------------------------- Adding Cyclotron ---------------------------- *" << endl;
gsell's avatar
gsell committed
354

355 356
    Cyclotron *elptr = dynamic_cast<Cyclotron *>(cycl.clone());
    myElements.push_back(elptr);
357

358
    // Is this a Spiral Inflector Simulation? If yes, we'll give the user some
359 360 361 362 363 364
    // useful information
    spiral_flag = elptr->getSpiralFlag();

    if(spiral_flag) {

        *gmsg << endl << "* This is a Spiral Inflector Simulation! This means the following:" << endl;
365
        *gmsg         << "* 1.) It is up to the user to provide appropriate geometry, electric and magnetic fields!" << endl;
366 367
        *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
368
        *gmsg         << "* 2.) For high currents it is strongly recommended to use the SAAMG fieldsolver," << endl;
369 370
        *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;
371 372
        *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;
373 374 375
        *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;
376
        mbHandler_m->setNumBunch(1);
377 378
    }

379
    // Fresh run (no restart):
380
    if(!OpalData::getInstance()->inRestartRun()) {
381

382
        // Get reference values from cyclotron element
383
        // For now, these are still stored in mm. should be the only ones. -DW
384 385 386 387
        referenceR     = elptr->getRinit();
        referenceTheta = elptr->getPHIinit();
        referenceZ     = elptr->getZinit();
        referencePr    = elptr->getPRinit();
388 389
        referencePz    = elptr->getPZinit();

390
        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
391 392
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                "PHIINIT is out of [-180, 180)!");
393 394 395 396 397
        }

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

        // Calculate reference azimuthal (tangential) momentum from total-, z- and radial momentum:
398
        float insqrt = referencePtot * referencePtot - \
399
            referencePr * referencePr - referencePz * referencePz;
400 401 402 403 404

        if(insqrt < 0) {

            if(insqrt > -1.0e-10) {

405
                referencePt = 0.0;
406 407 408

            } else {

409
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
410
                                    "Pt imaginary!");
411 412 413 414 415 416 417
            }

        } else {

            referencePt = sqrt(insqrt);
        }

418
        if(referencePtot < 0.0)
Daniel Winklehner's avatar
Daniel Winklehner committed
419
            referencePt *= -1.0;
420 421
        // End calculate referencePt

422
        // Restart a run:
423 424
    } else {

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

441
        // Adjust some of the reference variables from the h5 file
442 443
        referencePhi *= Physics::deg2rad;
        referencePsi *= Physics::deg2rad;
444
        referencePtot = bega;
445
        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
446
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
447 448
                                "PHIINIT is out of [-180, 180)!");
        }
449 450
    }

451 452
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
453

454
    *gmsg << endl;
adelmann's avatar
adelmann committed
455
    *gmsg << "* Bunch global starting position:" << endl;
456 457
    *gmsg << "* RINIT = " << referenceR  << " [mm]" << endl;
    *gmsg << "* PHIINIT = " << referenceTheta << " [deg]" << endl;
458
    *gmsg << "* ZINIT = " << referenceZ << " [mm]" << endl;
459
    *gmsg << endl;
adelmann's avatar
adelmann committed
460
    *gmsg << "* Bunch global starting momenta:" << endl;
gsell's avatar
gsell committed
461 462
    *gmsg << "* Initial gamma = " << itsReference.getGamma() << endl;
    *gmsg << "* Initial beta = " << itsReference.getBeta() << endl;
463
    *gmsg << "* Reference total momentum (beta * gamma) = " << referencePtot * 1000.0 << " [MCU]" << endl;
464 465 466
    *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;
467
    *gmsg << endl;
adelmann's avatar
adelmann committed
468

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

472 473 474
    // 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
475

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

479
    std::string type = elptr->getCyclotronType();
480
    *gmsg << "* Type of cyclotron = " << type << " " << endl;
481

482 483
    double rmin = elptr->getMinR();
    double rmax = elptr->getMaxR();
adelmann's avatar
adelmann committed
484
    *gmsg << "* Radial aperture = " << rmin << " ... " << rmax<<" [m] "<< endl;
485 486 487

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

    double h = elptr->getCyclHarm();
491
    *gmsg << "* Number of trimcoils = " << elptr->getNumberOfTrimcoils() << endl;
492
    *gmsg << "* Harmonic number h = " << h << " " << endl;
493

gsell's avatar
gsell committed
494 495 496 497 498 499 500
    /**
     * 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
501
     * fieldflag = 5, readin FFA format file for MSU/FNAL FFA
gsell's avatar
gsell committed
502
     * fieldflag = 6, readin both median plane B field map and 3D E field map of RF cavity for compact cyclotron
503
     * fieldflag = 7, read in fields for Daniel's synchrocyclotron simulations
gsell's avatar
gsell committed
504
     */
505
    int fieldflag = elptr->getFieldFlag(type);
gsell's avatar
gsell committed
506

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

snuverink_j's avatar
snuverink_j committed
510
    double BcParameter[8] = {};
511

512 513
    BcParameter[0] = 0.001 * elptr->getRmin();
    BcParameter[1] = 0.001 * elptr->getRmax();
gsell's avatar
gsell committed
514

515
    // Store inner radius and outer radius of cyclotron field map in the list
516
    buildupFieldList(BcParameter, ElementBase::CYCLOTRON, elptr);
gsell's avatar
gsell committed
517 518 519 520 521 522 523 524 525 526 527 528 529 530 531
}

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

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

534
    *gmsg << "* --------- Collimator -----------------------------" << endl;
gsell's avatar
gsell committed
535

536
    CCollimator* elptr = dynamic_cast<CCollimator *>(coll.clone());
537
    myElements.push_back(elptr);
gsell's avatar
gsell committed
538

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

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

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

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

551
    double zstart = elptr->getZStart();
adelmann's avatar
adelmann committed
552
    *gmsg << "* Zstart= " << zstart << " [mm]" << endl;
553 554

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

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

frey_m's avatar
frey_m committed
560
    elptr->initialise(itsBunch_m);
gsell's avatar
gsell committed
561

snuverink_j's avatar
snuverink_j committed
562
    double BcParameter[8] = {};
563

564 565 566 567 568 569
    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 ;

570
    buildupFieldList(BcParameter, ElementBase::CCOLLIMATOR, elptr);
gsell's avatar
gsell committed
571 572 573 574 575 576 577 578 579 580 581 582
}

/**
 *
 *
 * @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
583 584 585 586 587 588 589 590 591 592 593 594
/**
 *
 *
 * @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
595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614
/**
 *
 *
 * @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()));
}

615 616
/**
 *
617 618
 *
 *  @param
619 620 621 622 623
 */
void ParallelCyclotronTracker::visitFlexibleCollimator(const FlexibleCollimator &) {

}

gsell's avatar
gsell committed
624 625 626 627 628 629 630 631 632 633
/**
 *
 *
 * @param lamb
 */
void ParallelCyclotronTracker::visitLambertson(const Lambertson &lamb) {
    *gmsg << "In Lambertson; L= " << lamb.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Lambertson *>(lamb.clone()));
}

634 635 636
void ParallelCyclotronTracker::visitOffset(const Offset & off) {
    if (opalRing_m == NULL)
        throw OpalException(
637 638
                            "ParallelCylcotronTracker::visitOffset",
                            "Attempt to place an offset when Ring not defined");
639 640
    Offset* offNonConst = const_cast<Offset*>(&off);
    offNonConst->updateGeometry(opalRing_m->getNextPosition(),
641
                                opalRing_m->getNextNormal());
642 643 644
    opalRing_m->appendElement(off);
}

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

693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740
/**
 *
 *
 * @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
741 742 743 744 745 746
/**
 *
 *
 * @param prob
 */
void ParallelCyclotronTracker::visitProbe(const Probe &prob) {
747
    *gmsg << "* -----------  Probe -------------------------------" << endl;
748 749
    Probe *elptr = dynamic_cast<Probe *>(prob.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
750

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

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

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

760
    double yend = elptr->getYEnd();
761
    *gmsg << "YEnd= " << yend << " [mm]" << endl;
gsell's avatar
gsell committed
762 763

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

snuverink_j's avatar
snuverink_j committed
766
    double BcParameter[8] = {};
767

768 769 770 771
    BcParameter[0] = 0.001 * xstart ;
    BcParameter[1] = 0.001 * xend;
    BcParameter[2] = 0.001 * ystart ;
    BcParameter[3] = 0.001 * yend;
772
    BcParameter[4] = 0.001 ; // width
gsell's avatar
gsell committed
773 774

    // store probe parameters in the list
775
    buildupFieldList(BcParameter, ElementBase::PROBE, elptr);
gsell's avatar
gsell committed
776 777 778 779 780 781 782 783 784 785 786 787
}

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

788 789 790 791 792 793
void ParallelCyclotronTracker::visitSBend3D(const SBend3D &bend) {
    *gmsg << "Adding SBend3D" << endl;
    if (opalRing_m != NULL)
        opalRing_m->appendElement(bend);
    else
        throw OpalException("ParallelCyclotronTracker::visitSBend3D",
794
                            "Need to define a RINGDEFINITION to use SBend3D element");
795 796
}

ext-rogers_c's avatar
ext-rogers_c committed
797 798
void ParallelCyclotronTracker::visitScalingFFAMagnet(const ScalingFFAMagnet &bend) {
    *gmsg << "Adding ScalingFFAMagnet" << endl;
ext-rogers_c's avatar
ext-rogers_c committed
799
    if (opalRing_m != NULL) {
800
        opalRing_m->appendElement(bend);
ext-rogers_c's avatar
ext-rogers_c committed
801
    } else {
ext-rogers_c's avatar
ext-rogers_c committed
802 803
        throw OpalException("ParallelCyclotronTracker::visitScalingFFAMagnet",
                            "Need to define a RINGDEFINITION to use ScalingFFAMagnet element");
ext-rogers_c's avatar
ext-rogers_c committed
804
    }
805 806
}

807 808 809 810 811 812
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",
813
                            "Need to define a RINGDEFINITION to use VariableRFCavity element");
814 815 816 817 818 819 820 821 822 823
}

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

gsell's avatar
gsell committed
826 827 828 829 830 831 832 833
/**
 *
 *
 * @param as
 */
void ParallelCyclotronTracker::visitRFCavity(const RFCavity &as) {

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

835 836
    RFCavity *elptr = dynamic_cast<RFCavity *>(as.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
837 838 839 840 841 842 843

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

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

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

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

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

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

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

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

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

868 869 870
    /*
      Setup time dependence and in case of no
      timedependence use a polynom with  a_0 = 1 and a_k = 0, k = 1,2,3.
871
    */
872

873 874 875 876
    std::shared_ptr<AbstractTimeDependence> freq_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> ampl_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> phase_atd = nullptr;

877
    dvector_t  unityVec;
878 879 880 881
    unityVec.push_back(1.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
882

883
    if (elptr->getFrequencyModelName() != "") {
884 885
        freq_atd = AbstractTimeDependence::getTimeDependence(elptr->getFrequencyModelName());
        *gmsg << "* Variable frequency RF Model name " << elptr->getFrequencyModelName() << endl;
886
    }
887 888
    else
        freq_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
889 890

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

    if (elptr->getPhaseModelName() != "") {
898 899
        phase_atd = AbstractTimeDependence::getTimeDependence(elptr->getPhaseModelName());
        *gmsg << "* Variable phase RF Model name " << elptr->getPhaseModelName() << endl;
900
    }
901 902
    else
        phase_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
903

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

snuverink_j's avatar
snuverink_j committed
907
    double BcParameter[8] = {};
908

909 910
    BcParameter[0] = 0.001 * rmin;
    BcParameter[1] = 0.001 * rmax;
911
    BcParameter[2] = 0.001 * pdis;
gsell's avatar
gsell committed
912 913
    BcParameter[3] = angle;

914
    buildupFieldList(BcParameter, ElementBase::RFCAVITY, elptr);
gsell's avatar
gsell committed
915 916 917 918 919 920 921 922
}

/**
 *
 *
 * @param rfq
 */
void ParallelCyclotronTracker::visitRFQuadrupole(const RFQuadrupole &rfq) {
923
    *gmsg << "In RFQuadrupole; L = " << rfq.getElementLength() << " however the element is missing " << endl;
gsell's avatar
gsell committed
924 925 926 927 928 929 930 931 932
    myElements.push_back(dynamic_cast<RFQuadrupole *>(rfq.clone()));
}

/**
 *
 *
 * @param bend
 */
void ParallelCyclotronTracker::visitSBend(const SBend &bend) {
933
    *gmsg << "In SBend; L = " << bend.getElementLength() << " however the element is missing " << endl;
gsell's avatar
gsell committed
934 935 936 937 938 939 940 941 942 943 944 <