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

#include "Algorithms/ParallelCyclotronTracker.h"
20

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

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

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

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

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

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

#include "Physics/Physics.h"

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

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

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

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

86

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

89 90 91
Vector_t const ParallelCyclotronTracker::xaxis = Vector_t(1.0, 0.0, 0.0);
Vector_t const ParallelCyclotronTracker::yaxis = Vector_t(0.0, 1.0, 0.0);
Vector_t const ParallelCyclotronTracker::zaxis = Vector_t(0.0, 0.0, 1.0);
gsell's avatar
gsell committed
92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107

extern Inform *gmsg;

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

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

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

145 146 147 148 149 150 151 152 153
    // FIXME Change track command
    if ( initialTotalNum_m == 1 ) {
        mode_m = MODE::SINGLE;
    } else if ( initialTotalNum_m == 2 ) {
        mode_m = MODE::SEO;
    } else if ( initialTotalNum_m > 2 ) {
        mode_m = MODE::BUNCH;
    } else
        mode_m = MODE::UNDEFINED;
154

155 156 157 158 159 160 161 162
    if ( timeIntegrator == 0 ) {
        stepper_m = stepper::INTEGRATOR::RK4;
    } else if ( timeIntegrator == 1) {
        stepper_m = stepper::INTEGRATOR::LF2;
    } else if ( timeIntegrator == 2) {
        stepper_m = stepper::INTEGRATOR::MTS;
    } else
        stepper_m = stepper::INTEGRATOR::UNDEFINED;
gsell's avatar
gsell committed
163 164 165 166 167 168 169
}

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

frey_m's avatar
frey_m committed
180

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

189
    Inform msg("bgf_main_collision_test ");
190

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

195
    Vector_t intecoords = 0.0;
196

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

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

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


frey_m's avatar
frey_m committed
235 236 237 238 239 240
double ParallelCyclotronTracker::computeRadius(const Vector_t& meanR) const {
    // New OPAL 2.0: m --> mm
    return 1000.0 * std::sqrt(meanR(0) * meanR(0) + meanR(1) * meanR(1));
}


241 242 243 244
double ParallelCyclotronTracker::computePathLengthUpdate(const double& dt,
                                                         short bunchNr)
{
    double dotP = 0.0;
245
    if ( Options::psDumpFrame == Options::BUNCH_MEAN || isMultiBunch()) {
246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269
        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
}


270

gsell's avatar
gsell committed
271 272 273 274 275
/**
 *
 *
 * @param fn Base file name
 */
276
void ParallelCyclotronTracker::openFiles(size_t numFiles, std::string SfileName) {
gsell's avatar
gsell committed
277

278
    for (unsigned int i=0; i<numFiles; i++) {
279 280 281 282 283 284 285 286
        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
287

288 289 290 291 292 293
        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;
294
    }
gsell's avatar
gsell committed
295 296 297 298 299 300 301 302
}

/**
 * Close all files related to
 * special output in the Cyclotron
 * mode.
 */
void ParallelCyclotronTracker::closeFiles() {
303
    for (auto & file : outfTheta_m) {
304
        file->close();
305
    }
gsell's avatar
gsell committed
306 307
}

308
/**
309 310 311
 *
 * @param ring
 */
kraus's avatar
kraus committed
312
void ParallelCyclotronTracker::visitRing(const Ring &ring) {
Daniel Winklehner's avatar
Daniel Winklehner committed
313

314
    *gmsg << "* ----------------------------- Adding Ring ------------------------------ *" << endl;
Daniel Winklehner's avatar
Daniel Winklehner committed
315

snuverink_j's avatar
snuverink_j committed
316
    delete opalRing_m;
Daniel Winklehner's avatar
Daniel Winklehner committed
317

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

320
    myElements.push_back(opalRing_m);
Daniel Winklehner's avatar
Daniel Winklehner committed
321

frey_m's avatar
frey_m committed
322
    opalRing_m->initialise(itsBunch_m);
323 324 325 326

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

328
    if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
kraus's avatar
kraus committed
329
        throw OpalException("Error in ParallelCyclotronTracker::visitRing",
330 331
                            "PHIINIT is out of [-180, 180)!");
    }
Daniel Winklehner's avatar
Daniel Winklehner committed
332 333

    referenceZ = 0.0;
334
    referencePz = 0.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
335 336 337 338

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

339 340
    if(referencePtot < 0.0)
        referencePt *= -1.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
341

342 343
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
gsell's avatar
gsell committed
344

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

347
    buildupFieldList(BcParameter, ElementBase::RING, opalRing_m);
348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363

    // 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
364 365 366 367 368 369 370 371

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

372
    *gmsg << "* -------------------------- Adding Cyclotron ---------------------------- *" << endl;
gsell's avatar
gsell committed
373

374 375
    Cyclotron *elptr = dynamic_cast<Cyclotron *>(cycl.clone());
    myElements.push_back(elptr);
376

377
    // Is this a Spiral Inflector Simulation? If yes, we'll give the user some
378 379 380 381 382 383
    // useful information
    spiral_flag = elptr->getSpiralFlag();

    if(spiral_flag) {

        *gmsg << endl << "* This is a Spiral Inflector Simulation! This means the following:" << endl;
384
        *gmsg         << "* 1.) It is up to the user to provide appropriate geometry, electric and magnetic fields!" << endl;
385 386
        *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
387
        *gmsg         << "* 2.) For high currents it is strongly recommended to use the SAAMG fieldsolver," << endl;
388 389
        *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;
390 391
        *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;
392 393 394
        *gmsg         << "*     (This is not an issue for spiral inflectors as they are typically < 100 keV/amu.)" << endl;
        *gmsg << endl << "* Note: For now, multi-bunch mode (MBM) needs to be de-activated for spiral inflector" << endl;
        *gmsg         << "* and space charge needs to be solved every time-step. numBunch_m and scSolveFreq are reset." << endl;
frey_m's avatar
frey_m committed
395 396 397
        if (isMultiBunch()) {
            mbHandler_m = nullptr;
        }
398 399
    }

400
    // Fresh run (no restart):
401
    if(!OpalData::getInstance()->inRestartRun()) {
402

403
        // Get reference values from cyclotron element
404
        // For now, these are still stored in mm. should be the only ones. -DW
405 406 407 408
        referenceR     = elptr->getRinit();
        referenceTheta = elptr->getPHIinit();
        referenceZ     = elptr->getZinit();
        referencePr    = elptr->getPRinit();
409 410
        referencePz    = elptr->getPZinit();

411
        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
412 413
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                "PHIINIT is out of [-180, 180)!");
414 415 416 417 418
        }

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

        // Calculate reference azimuthal (tangential) momentum from total-, z- and radial momentum:
419
        float insqrt = referencePtot * referencePtot - \
420
            referencePr * referencePr - referencePz * referencePz;
421 422 423 424 425

        if(insqrt < 0) {

            if(insqrt > -1.0e-10) {

426
                referencePt = 0.0;
427 428 429

            } else {

430
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
431
                                    "Pt imaginary!");
432 433 434 435 436 437 438
            }

        } else {

            referencePt = sqrt(insqrt);
        }

439
        if(referencePtot < 0.0)
Daniel Winklehner's avatar
Daniel Winklehner committed
440
            referencePt *= -1.0;
441 442
        // End calculate referencePt

443
        // Restart a run:
444 445
    } else {

446
        // If the user wants to save the restarted run in local frame,
447
        // make sure the previous h5 file was local too
448
      if (Options::psDumpFrame != Options::GLOBAL) {
winklehner_d's avatar
-DW  
winklehner_d committed
449
        if (!previousH5Local) {
450 451
                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
452
        }
453 454
            // 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
455 456
      } else {
        if (previousH5Local) {
457 458
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                    "You are trying a global restart from a local h5 file!");
459
            }
460
        }
461

462
        // Adjust some of the reference variables from the h5 file
463 464
        referencePhi *= Physics::deg2rad;
        referencePsi *= Physics::deg2rad;
465
        referencePtot = bega;
466
        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
467
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
468 469
                                "PHIINIT is out of [-180, 180)!");
        }
470 471
    }

472 473
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
474

475
    *gmsg << endl;
adelmann's avatar
adelmann committed
476
    *gmsg << "* Bunch global starting position:" << endl;
477 478
    *gmsg << "* RINIT = " << referenceR  << " [mm]" << endl;
    *gmsg << "* PHIINIT = " << referenceTheta << " [deg]" << endl;
479
    *gmsg << "* ZINIT = " << referenceZ << " [mm]" << endl;
480
    *gmsg << endl;
adelmann's avatar
adelmann committed
481
    *gmsg << "* Bunch global starting momenta:" << endl;
gsell's avatar
gsell committed
482 483
    *gmsg << "* Initial gamma = " << itsReference.getGamma() << endl;
    *gmsg << "* Initial beta = " << itsReference.getBeta() << endl;
484
    *gmsg << "* Reference total momentum (beta * gamma) = " << referencePtot * 1000.0 << " [MCU]" << endl;
485 486 487
    *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;
488
    *gmsg << endl;
adelmann's avatar
adelmann committed
489

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

493 494 495
    // 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
496

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

500
    std::string type = elptr->getCyclotronType();
501
    *gmsg << "* Type of cyclotron = " << type << " " << endl;
502

503 504
    double rmin = elptr->getMinR();
    double rmax = elptr->getMaxR();
adelmann's avatar
adelmann committed
505
    *gmsg << "* Radial aperture = " << rmin << " ... " << rmax<<" [m] "<< endl;
506 507 508

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

    double h = elptr->getCyclHarm();
512
    *gmsg << "* Number of trimcoils = " << elptr->getNumberOfTrimcoils() << endl;
513
    *gmsg << "* Harmonic number h = " << h << " " << endl;
514

gsell's avatar
gsell committed
515 516 517 518 519 520 521
    /**
     * 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
522
     * fieldflag = 5, readin FFA format file for MSU/FNAL FFA
gsell's avatar
gsell committed
523
     * fieldflag = 6, readin both median plane B field map and 3D E field map of RF cavity for compact cyclotron
524
     * fieldflag = 7, read in fields for Daniel's synchrocyclotron simulations
gsell's avatar
gsell committed
525
     */
526
    int fieldflag = elptr->getFieldFlag(type);
gsell's avatar
gsell committed
527

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

snuverink_j's avatar
snuverink_j committed
531
    double BcParameter[8] = {};
532

533 534
    BcParameter[0] = 0.001 * elptr->getRmin();
    BcParameter[1] = 0.001 * elptr->getRmax();
gsell's avatar
gsell committed
535

536
    // Store inner radius and outer radius of cyclotron field map in the list
537
    buildupFieldList(BcParameter, ElementBase::CYCLOTRON, elptr);
gsell's avatar
gsell committed
538 539 540 541 542 543 544 545 546 547 548 549 550 551 552
}

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

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

555
    *gmsg << "* --------- Collimator -----------------------------" << endl;
gsell's avatar
gsell committed
556

557
    CCollimator* elptr = dynamic_cast<CCollimator *>(coll.clone());
558
    myElements.push_back(elptr);
gsell's avatar
gsell committed
559

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

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

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

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

572
    double zstart = elptr->getZStart();
adelmann's avatar
adelmann committed
573
    *gmsg << "* Zstart= " << zstart << " [mm]" << endl;
574 575

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

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

frey_m's avatar
frey_m committed
581
    elptr->initialise(itsBunch_m);
gsell's avatar
gsell committed
582

snuverink_j's avatar
snuverink_j committed
583
    double BcParameter[8] = {};
584

585 586 587 588 589 590
    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 ;

591
    buildupFieldList(BcParameter, ElementBase::CCOLLIMATOR, elptr);
gsell's avatar
gsell committed
592 593 594 595 596 597 598 599 600 601 602 603
}

/**
 *
 *
 * @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
604 605 606 607 608 609 610 611 612 613 614 615
/**
 *
 *
 * @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
616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635
/**
 *
 *
 * @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()));
}

636 637
/**
 *
638 639
 *
 *  @param
640 641 642 643 644
 */
void ParallelCyclotronTracker::visitFlexibleCollimator(const FlexibleCollimator &) {

}

gsell's avatar
gsell committed
645 646 647 648 649 650 651 652 653 654
/**
 *
 *
 * @param lamb
 */
void ParallelCyclotronTracker::visitLambertson(const Lambertson &lamb) {
    *gmsg << "In Lambertson; L= " << lamb.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Lambertson *>(lamb.clone()));
}

655 656 657
void ParallelCyclotronTracker::visitOffset(const Offset & off) {
    if (opalRing_m == NULL)
        throw OpalException(
658 659
                            "ParallelCylcotronTracker::visitOffset",
                            "Attempt to place an offset when Ring not defined");
660 661
    Offset* offNonConst = const_cast<Offset*>(&off);
    offNonConst->updateGeometry(opalRing_m->getNextPosition(),
662
                                opalRing_m->getNextNormal());
663 664 665
    opalRing_m->appendElement(off);
}

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

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 749 750 751 752 753 754 755 756 757 758 759 760 761
/**
 *
 *
 * @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
762 763 764 765 766 767
/**
 *
 *
 * @param prob
 */
void ParallelCyclotronTracker::visitProbe(const Probe &prob) {
768
    *gmsg << "* -----------  Probe -------------------------------" << endl;
769 770
    Probe *elptr = dynamic_cast<Probe *>(prob.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
771

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

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

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

781
    double yend = elptr->getYEnd();
782
    *gmsg << "YEnd= " << yend << " [mm]" << endl;
gsell's avatar
gsell committed
783 784

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

snuverink_j's avatar
snuverink_j committed
787
    double BcParameter[8] = {};
788

789 790 791 792
    BcParameter[0] = 0.001 * xstart ;
    BcParameter[1] = 0.001 * xend;
    BcParameter[2] = 0.001 * ystart ;
    BcParameter[3] = 0.001 * yend;
793
    BcParameter[4] = 0.001 ; // width
gsell's avatar
gsell committed
794 795

    // store probe parameters in the list
796
    buildupFieldList(BcParameter, ElementBase::PROBE, elptr);
gsell's avatar
gsell committed
797 798 799 800 801 802 803 804 805 806 807 808
}

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

809 810 811 812 813 814
void ParallelCyclotronTracker::visitSBend3D(const SBend3D &bend) {
    *gmsg << "Adding SBend3D" << endl;
    if (opalRing_m != NULL)
        opalRing_m->appendElement(bend);
    else
        throw OpalException("ParallelCyclotronTracker::visitSBend3D",
815
                            "Need to define a RINGDEFINITION to use SBend3D element");
816 817
}

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

828 829 830 831 832 833
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",
834
                            "Need to define a RINGDEFINITION to use VariableRFCavity element");
835 836 837 838 839 840 841 842 843 844
}

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

gsell's avatar
gsell committed
847 848 849 850 851 852 853 854
/**
 *
 *
 * @param as
 */
void ParallelCyclotronTracker::visitRFCavity(const RFCavity &as) {

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

856 857
    RFCavity *elptr = dynamic_cast<RFCavity *>(as.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
858 859 860 861 862 863 864

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

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

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

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

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

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

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

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

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

889 890 891
    /*
      Setup time dependence and in case of no
      timedependence use a polynom with  a_0 = 1 and a_k = 0, k = 1,2,3.
892
    */
893

894 895 896 897
    std::shared_ptr<AbstractTimeDependence> freq_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> ampl_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> phase_atd = nullptr;

898
    dvector_t  unityVec;
899 900 901 902
    unityVec.push_back(1.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
903

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

    if (elptr->getAmplitudeModelName() != "") {
912 913
        ampl_atd = AbstractTimeDependence::getTimeDependence(elptr->getAmplitudeModelName());
        *gmsg << "* Variable amplitude RF Model name " << elptr->getAmplitudeModelName() << endl;
914
    }
915 916
    else
        ampl_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
917 918

    if (elptr->getPhaseModelName() != "") {
919 920
        phase_atd = AbstractTimeDependence::getTimeDependence(elptr->getPhaseModelName());
        *gmsg << "* Variable phase RF Model name " << elptr->getPhaseModelName() << endl;
921
    }
922 923
    else
        phase_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
924

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