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 120 121 122 123 124 125 126
    : Tracker(beamline, bunch, reference, revBeam, revTrack)
    , bgf_m(nullptr)
    , maxSteps_m(maxSTEPS)
    , lastDumpedStep_m(0)
    , myNode_m(Ippl::myNode())
    , initialLocalNum_m(bunch->getLocalNum())
    , initialTotalNum_m(bunch->getTotalNum())
    , opalRing_m(nullptr)
    , itsStepper_mp(nullptr)
{
gsell's avatar
gsell committed
127 128 129
    itsBeamline = dynamic_cast<Beamline *>(beamline.clone());
    itsDataSink = &ds;

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

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

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

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

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

frey_m's avatar
frey_m committed
179

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

188
    Inform msg("bgf_main_collision_test ");
189

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

194
    Vector_t intecoords = 0.0;
195

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

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

209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226
// 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;
}


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


256

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

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

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

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

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

300
    *gmsg << "* ----------------------------- Adding Ring ------------------------------ *" << endl;
Daniel Winklehner's avatar
Daniel Winklehner committed
301

snuverink_j's avatar
snuverink_j committed
302
    delete opalRing_m;
Daniel Winklehner's avatar
Daniel Winklehner committed
303

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

306
    myElements.push_back(opalRing_m);
Daniel Winklehner's avatar
Daniel Winklehner committed
307

frey_m's avatar
frey_m committed
308
    opalRing_m->initialise(itsBunch_m);
309 310 311 312

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

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

    referenceZ = 0.0;
320
    referencePz = 0.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
321 322 323 324

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

325 326
    if(referencePtot < 0.0)
        referencePt *= -1.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
327

328 329
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
gsell's avatar
gsell committed
330

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

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

    // 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
350 351 352 353 354 355 356 357

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

358
    *gmsg << "* -------------------------- Adding Cyclotron ---------------------------- *" << endl;
gsell's avatar
gsell committed
359

360 361
    Cyclotron *elptr = dynamic_cast<Cyclotron *>(cycl.clone());
    myElements.push_back(elptr);
362

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

    if(spiral_flag) {

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

384
    // Fresh run (no restart):
385
    if(!OpalData::getInstance()->inRestartRun()) {
386

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

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

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

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

        if(insqrt < 0) {

            if(insqrt > -1.0e-10) {

410
                referencePt = 0.0;
411 412 413

            } else {

414
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
415
                                    "Pt imaginary!");
416 417 418 419 420 421 422
            }

        } else {

            referencePt = sqrt(insqrt);
        }

423
        if(referencePtot < 0.0)
Daniel Winklehner's avatar
Daniel Winklehner committed
424
            referencePt *= -1.0;
425 426
        // End calculate referencePt

427
        // Restart a run:
428 429
    } else {

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

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

456 457
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
458

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

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

477 478 479
    // 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
480

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

484
    std::string type = elptr->getCyclotronType();
485
    *gmsg << "* Type of cyclotron = " << type << " " << endl;
486

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

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

    double h = elptr->getCyclHarm();
496
    *gmsg << "* Number of trimcoils = " << elptr->getNumberOfTrimcoils() << endl;
497
    *gmsg << "* Harmonic number h = " << h << " " << endl;
498

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

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

snuverink_j's avatar
snuverink_j committed
515
    double BcParameter[8] = {};
516

517 518
    BcParameter[0] = 0.001 * elptr->getRmin();
    BcParameter[1] = 0.001 * elptr->getRmax();
gsell's avatar
gsell committed
519

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

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

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

539
    *gmsg << "* --------- Collimator -----------------------------" << endl;
gsell's avatar
gsell committed
540

541
    CCollimator* elptr = dynamic_cast<CCollimator *>(coll.clone());
542
    myElements.push_back(elptr);
gsell's avatar
gsell committed
543

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

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

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

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

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

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

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

frey_m's avatar
frey_m committed
565
    elptr->initialise(itsBunch_m);
gsell's avatar
gsell committed
566

snuverink_j's avatar
snuverink_j committed
567
    double BcParameter[8] = {};
568

569 570 571 572 573 574
    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 ;

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

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

620 621
/**
 *
622 623
 *
 *  @param
624 625 626 627 628
 */
void ParallelCyclotronTracker::visitFlexibleCollimator(const FlexibleCollimator &) {

}

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

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

gsell's avatar
gsell committed
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 677 678 679 680 681
/**
 *
 *
 * @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
682 683 684 685 686 687 688
/**
 *
 *
 * @param multT
 */
void ParallelCyclotronTracker::visitMultipoleT(const MultipoleT &multT) {
    *gmsg << "Adding MultipoleT" << endl;
689
    if (opalRing_m != NULL) {
ext-rogers_c's avatar
ext-rogers_c committed
690
        opalRing_m->appendElement(multT);
691
    } else {
ext-rogers_c's avatar
ext-rogers_c committed
692 693
        throw OpalException("ParallelCyclotronTracker::visitMultipoleT",
                            "Need to define a RINGDEFINITION to use MultipoleT element");
694
    }
ext-rogers_c's avatar
ext-rogers_c committed
695 696 697
    myElements.push_back(dynamic_cast<MultipoleT *>(multT.clone()));
}

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 741 742 743 744 745
/**
 *
 *
 * @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
746 747 748 749 750 751
/**
 *
 *
 * @param prob
 */
void ParallelCyclotronTracker::visitProbe(const Probe &prob) {
752
    *gmsg << "* -----------  Probe -------------------------------" << endl;
753 754
    Probe *elptr = dynamic_cast<Probe *>(prob.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
755

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

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

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

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

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

snuverink_j's avatar
snuverink_j committed
771
    double BcParameter[8] = {};
772

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

878 879 880 881
    std::shared_ptr<AbstractTimeDependence> freq_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> ampl_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> phase_atd = nullptr;

882
    dvector_t  unityVec;
883 884 885 886
    unityVec.push_back(1.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
887

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

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

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

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

snuverink_j's avatar
snuverink_j committed
912
    double BcParameter[8] = {};
913

914 915
    BcParameter[0] = 0.001 * rmin;
    BcParameter[1] = 0.001 * rmax;
916
    BcParameter[2] = 0.001 * pdis;
gsell's avatar
gsell committed
917 918
    BcParameter[3] = angle;

919
    buildupFieldList(BcParameter, ElementBase::RFCAVITY, elptr);
gsell's avatar
gsell committed
920 921 922 923 924 925 926 927
}

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