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

#include "Algorithms/ParallelCyclotronTracker.h"
20

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

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

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

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

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

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

#include "Physics/Physics.h"

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

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

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

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

86

87
constexpr 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
void ParallelCyclotronTracker::bgf_main_collision_test() {
182
    if(!bgf_m) return;
183

184
    Inform msg("bgf_main_collision_test ");
185

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

190
    Vector_t intecoords = 0.0;
191

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

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

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


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


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

frey_m's avatar
frey_m committed
244 245
        for(unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) {
            dotP[itsBunch_m->bunchNum[i]] += dot(itsBunch_m->P[i], itsBunch_m->P[i]);
246 247
        }

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

frey_m's avatar
frey_m committed
250 251 252
        // dot-product over all particles
        double sum = std::accumulate(dotP.begin(), dotP.end() - 1, 0);
        dotP.back() = sum / double(itsBunch_m->getTotalNum());
frey_m's avatar
frey_m committed
253

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

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

frey_m's avatar
frey_m committed
267 268 269 270 271
    for (size_t i = 0; i < dotP.size(); ++i) {
        double const gamma = std::sqrt(1.0 + dotP[i]);
        double const beta  = std::sqrt(dotP[i]) / gamma;
        dl[i] = c_mmtns * dt * 1.0e-3 * beta; // unit: m
    }
272 273 274
}


275

gsell's avatar
gsell committed
276 277 278 279 280
/**
 *
 *
 * @param fn Base file name
 */
281
void ParallelCyclotronTracker::openFiles(size_t numFiles, std::string SfileName) {
gsell's avatar
gsell committed
282

283
    for (unsigned int i=0; i<numFiles; i++) {
284 285 286 287 288 289 290 291
        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
292

293 294 295 296 297 298
        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;
299
    }
gsell's avatar
gsell committed
300 301 302 303 304 305 306 307
}

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

313
/**
314 315 316
 *
 * @param ring
 */
kraus's avatar
kraus committed
317
void ParallelCyclotronTracker::visitRing(const Ring &ring) {
Daniel Winklehner's avatar
Daniel Winklehner committed
318

319
    *gmsg << "* ----------------------------- Adding Ring ------------------------------ *" << endl;
Daniel Winklehner's avatar
Daniel Winklehner committed
320

snuverink_j's avatar
snuverink_j committed
321
    delete opalRing_m;
Daniel Winklehner's avatar
Daniel Winklehner committed
322

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

325
    myElements.push_back(opalRing_m);
Daniel Winklehner's avatar
Daniel Winklehner committed
326

frey_m's avatar
frey_m committed
327
    opalRing_m->initialise(itsBunch_m);
328 329 330 331

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

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

    referenceZ = 0.0;
339
    referencePz = 0.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
340 341 342 343

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

344 345
    if(referencePtot < 0.0)
        referencePt *= -1.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
346

347 348
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
gsell's avatar
gsell committed
349

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

352
    buildupFieldList(BcParameter, ElementBase::RING, opalRing_m);
353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368

    // 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
369 370 371 372 373 374 375 376

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

377
    *gmsg << "* -------------------------- Adding Cyclotron ---------------------------- *" << endl;
gsell's avatar
gsell committed
378

379 380
    Cyclotron *elptr = dynamic_cast<Cyclotron *>(cycl.clone());
    myElements.push_back(elptr);
381

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

    if(spiral_flag) {

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

405
    // Fresh run (no restart):
406
    if(!OpalData::getInstance()->inRestartRun()) {
407

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

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

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

        // Calculate reference azimuthal (tangential) momentum from total-, z- and radial momentum:
424
        float insqrt = referencePtot * referencePtot - \
425
            referencePr * referencePr - referencePz * referencePz;
426 427 428 429 430

        if(insqrt < 0) {

            if(insqrt > -1.0e-10) {

431
                referencePt = 0.0;
432 433 434

            } else {

435
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
436
                                    "Pt imaginary!");
437 438 439 440 441 442 443
            }

        } else {

            referencePt = sqrt(insqrt);
        }

444
        if(referencePtot < 0.0)
Daniel Winklehner's avatar
Daniel Winklehner committed
445
            referencePt *= -1.0;
446 447
        // End calculate referencePt

448
        // Restart a run:
449 450
    } else {

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

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

477 478
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
479

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

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

498 499 500
    // 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
501

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

505
    std::string type = elptr->getCyclotronType();
506
    *gmsg << "* Type of cyclotron = " << type << " " << endl;
507

508 509
    double rmin = elptr->getMinR();
    double rmax = elptr->getMaxR();
adelmann's avatar
adelmann committed
510
    *gmsg << "* Radial aperture = " << rmin << " ... " << rmax<<" [m] "<< endl;
511 512 513

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

    double h = elptr->getCyclHarm();
517
    *gmsg << "* Number of trimcoils = " << elptr->getNumberOfTrimcoils() << endl;
518
    *gmsg << "* Harmonic number h = " << h << " " << endl;
519

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

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

snuverink_j's avatar
snuverink_j committed
536
    double BcParameter[8] = {};
537

538 539
    BcParameter[0] = 0.001 * elptr->getRmin();
    BcParameter[1] = 0.001 * elptr->getRmax();
gsell's avatar
gsell committed
540

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

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

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

560
    *gmsg << "* --------- Collimator -----------------------------" << endl;
gsell's avatar
gsell committed
561

562
    CCollimator* elptr = dynamic_cast<CCollimator *>(coll.clone());
563
    myElements.push_back(elptr);
gsell's avatar
gsell committed
564

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

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

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

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

577
    double zstart = elptr->getZStart();
adelmann's avatar
adelmann committed
578
    *gmsg << "* Zstart= " << zstart << " [mm]" << endl;
579 580

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

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

frey_m's avatar
frey_m committed
586
    elptr->initialise(itsBunch_m);
gsell's avatar
gsell committed
587

snuverink_j's avatar
snuverink_j committed
588
    double BcParameter[8] = {};
589

590 591 592 593 594 595
    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 ;

596
    buildupFieldList(BcParameter, ElementBase::CCOLLIMATOR, elptr);
gsell's avatar
gsell committed
597 598 599 600 601 602 603 604 605 606 607 608
}

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

641 642
/**
 *
643 644
 *
 *  @param
645 646 647 648 649
 */
void ParallelCyclotronTracker::visitFlexibleCollimator(const FlexibleCollimator &) {

}

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

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

gsell's avatar
gsell committed
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 698 699 700 701 702
/**
 *
 *
 * @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
703 704 705 706 707 708 709
/**
 *
 *
 * @param multT
 */
void ParallelCyclotronTracker::visitMultipoleT(const MultipoleT &multT) {
    *gmsg << "Adding MultipoleT" << endl;
710
    if (opalRing_m != NULL) {
ext-rogers_c's avatar
ext-rogers_c committed
711
        opalRing_m->appendElement(multT);
712
    } else {
ext-rogers_c's avatar
ext-rogers_c committed
713 714
        throw OpalException("ParallelCyclotronTracker::visitMultipoleT",
                            "Need to define a RINGDEFINITION to use MultipoleT element");
715
    }
ext-rogers_c's avatar
ext-rogers_c committed
716 717 718
    myElements.push_back(dynamic_cast<MultipoleT *>(multT.clone()));
}

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 762 763 764 765 766
/**
 *
 *
 * @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
767 768 769 770 771 772
/**
 *
 *
 * @param prob
 */
void ParallelCyclotronTracker::visitProbe(const Probe &prob) {
773
    *gmsg << "* -----------  Probe -------------------------------" << endl;
774 775
    Probe *elptr = dynamic_cast<Probe *>(prob.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
776

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

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

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

786
    double yend = elptr->getYEnd();
787
    *gmsg << "YEnd= " << yend << " [mm]" << endl;
gsell's avatar
gsell committed
788 789

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

snuverink_j's avatar
snuverink_j committed
792
    double BcParameter[8] = {};
793

794 795 796 797
    BcParameter[0] = 0.001 * xstart ;
    BcParameter[1] = 0.001 * xend;
    BcParameter[2] = 0.001 * ystart ;
    BcParameter[3] = 0.001 * yend;
798
    BcParameter[4] = 0.001 ; // width
gsell's avatar
gsell committed
799 800

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

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

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

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

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

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

gsell's avatar
gsell committed
852 853 854 855 856 857 858 859
/**
 *
 *
 * @param as
 */
void ParallelCyclotronTracker::visitRFCavity(const RFCavity &as) {

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

861 862
    RFCavity *elptr = dynamic_cast<RFCavity *>(as.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
863 864 865 866 867 868 869

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

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

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

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

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

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

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

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

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

894 895 896
    /*
      Setup time dependence and in case of no
      timedependence use a polynom with  a_0 = 1 and a_k = 0, k = 1,2,3.
897
    */
898

899 900 901 902
    std::shared_ptr<AbstractTimeDependence> freq_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> ampl_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> phase_atd = nullptr;

903
    dvector_t  unityVec;
904 905 906 907
    unityVec.push_back(1.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);