ParallelCyclotronTracker.cpp 127 KB
Newer Older
gsell's avatar
gsell committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14
// ------------------------------------------------------------------------
// $RCSfile: ParallelCyclotronTracker.cpp,v $
// ------------------------------------------------------------------------
// $Revision: 1.1 $initialLocalNum_m
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: ParallelCyclotronTracker
//   The class for tracking particles with 3D space charge in Cyclotrons and FFAG's
//
// ------------------------------------------------------------------------
//
// $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/Ctunes.hpp>
kraus's avatar
kraus committed
20
#include "Algorithms/ParallelCyclotronTracker.h"
21 22
#include "Algorithms/PolynomialTimeDependence.h"
#include "Elements/OpalPolynomialTimeDependence.h"
23

gsell's avatar
gsell committed
24 25 26 27
#include <cfloat>
#include <iostream>
#include <fstream>
#include <vector>
28
#include "AbstractObjects/OpalData.h"
gsell's avatar
gsell committed
29 30 31 32

#include "AbsBeamline/Collimator.h"
#include "AbsBeamline/Corrector.h"
#include "AbsBeamline/Cyclotron.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"
gsell's avatar
gsell committed
43 44 45 46 47
#include "AbsBeamline/Probe.h"
#include "AbsBeamline/RBend.h"
#include "AbsBeamline/RFCavity.h"
#include "AbsBeamline/RFQuadrupole.h"
#include "AbsBeamline/SBend.h"
48
#include "AbsBeamline/SBend3D.h"
49
#include "AbsBeamline/ScalingFFAGMagnet.h"
gsell's avatar
gsell committed
50 51 52 53 54
#include "AbsBeamline/Separator.h"
#include "AbsBeamline/Septum.h"
#include "AbsBeamline/Solenoid.h"
#include "AbsBeamline/CyclotronValley.h"
#include "AbsBeamline/Stripper.h"
55
#include "AbsBeamline/VariableRFCavity.h"
56

57 58
#include "AbstractObjects/Element.h"

59
#include "Elements/OpalBeamline.h"
kraus's avatar
kraus committed
60
#include "AbsBeamline/Ring.h"
gsell's avatar
gsell committed
61 62 63

#include "BeamlineGeometry/Euclid3D.h"
#include "BeamlineGeometry/PlanarArcGeometry.h"
Jianjun Yang's avatar
Jianjun Yang committed
64
#include "BeamlineGeometry/RBendGeometry.h"
ext-rogers_c's avatar
ext-rogers_c committed
65
#include "BeamlineGeometry/StraightGeometry.h"
gsell's avatar
gsell committed
66 67 68 69 70 71 72 73 74 75 76 77
#include "Beamlines/Beamline.h"

#include "Fields/BMultipoleField.h"
#include "FixedAlgebra/FTps.h"
#include "FixedAlgebra/FTpsMath.h"
#include "FixedAlgebra/FVps.h"

#include "Physics/Physics.h"

#include "Utilities/NumToStr.h"
#include "Utilities/OpalException.h"

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

81
#include "Structure/H5PartWrapperForPC.h"
82
#include "Structure/BoundaryGeometry.h"
83
#include "Utilities/Options.h"
gsell's avatar
gsell committed
84 85 86 87 88 89 90

#include "Ctunes.h"
#include <cassert>

#include <hdf5.h>
#include "H5hut.h"

frey_m's avatar
frey_m committed
91 92
//FIXME Remove headers and dynamic_cast in readOneBunchFromFile
#include "Algorithms/PartBunch.h"
93
#ifdef ENABLE_AMR
frey_m's avatar
frey_m committed
94 95
    #include "Algorithms/AmrPartBunch.h"
#endif
frey_m's avatar
frey_m committed
96

gsell's avatar
gsell committed
97 98
class Beamline;
class PartData;
99

gsell's avatar
gsell committed
100 101 102
using Physics::pi;
using Physics::q_e;

103

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

106 107 108
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
109

snuverink_j's avatar
snuverink_j committed
110
//#define PSdim 6
gsell's avatar
gsell committed
111 112 113 114 115 116 117 118 119 120 121 122 123 124

extern Inform *gmsg;

// typedef FVector<double, PSdim> Vector;

/**
 * Constructor ParallelCyclotronTracker
 *
 * @param beamline
 * @param reference
 * @param revBeam
 * @param revTrack
 */
ParallelCyclotronTracker::ParallelCyclotronTracker(const Beamline &beamline,
125 126
                                                   const PartData &reference,
                                                   bool revBeam, bool revTrack):
gsell's avatar
gsell committed
127
    Tracker(beamline, reference, revBeam, revTrack),
128 129
    itsDataSink(nullptr),
    bgf_m(nullptr),
130
    lastDumpedStep_m(0),
131
    eta_m(0.01),
132 133
    initialR_m(nullptr),
    initialP_m(nullptr),
gsell's avatar
gsell committed
134 135
    myNode_m(Ippl::myNode()),
    initialLocalNum_m(0),
136
    initialTotalNum_m(0),
137
    opalRing_m(nullptr),
138 139
    itsStepper_mp(nullptr),
    mode_m(MODE::UNDEFINED),
140
    stepper_m(stepper::INTEGRATOR::UNDEFINED) {
gsell's avatar
gsell committed
141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156
    itsBeamline = dynamic_cast<Beamline *>(beamline.clone());
}

/**
 * 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
157
                                                   PartBunchBase<double, 3> *bunch,
gsell's avatar
gsell committed
158 159 160
                                                   DataSink &ds,
                                                   const PartData &reference,
                                                   bool revBeam, bool revTrack,
161
                                                   int maxSTEPS, int timeIntegrator):
162
    Tracker(beamline, bunch, reference, revBeam, revTrack),
163
    bgf_m(nullptr),
gsell's avatar
gsell committed
164
    maxSteps_m(maxSTEPS),
165
    lastDumpedStep_m(0),
166
    eta_m(0.01),
167 168
    initialR_m(nullptr),
    initialP_m(nullptr),
gsell's avatar
gsell committed
169
    myNode_m(Ippl::myNode()),
frey_m's avatar
frey_m committed
170 171
    initialLocalNum_m(bunch->getLocalNum()),
    initialTotalNum_m(bunch->getTotalNum()),
172 173
    opalRing_m(nullptr),
    itsStepper_mp(nullptr) {
gsell's avatar
gsell committed
174 175
    itsBeamline = dynamic_cast<Beamline *>(beamline.clone());
    itsDataSink = &ds;
frey_m's avatar
frey_m committed
176
    //  scaleFactor_m = itsBunch_m->getdT() * c;
gsell's avatar
gsell committed
177 178 179 180 181 182 183
    scaleFactor_m = 1;
    multiBunchMode_m = 0;

    IntegrationTimer_m = IpplTimings::getTimer("Integration");
    TransformTimer_m   = IpplTimings::getTimer("Frametransform");
    DumpTimer_m        = IpplTimings::getTimer("Dump");
    BinRepartTimer_m   = IpplTimings::getTimer("Binaryrepart");
184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202
    
    // 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;
    
    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
203 204 205 206 207 208 209
}

/**
 * Destructor ParallelCyclotronTracker
 *
 */
ParallelCyclotronTracker::~ParallelCyclotronTracker() {
210 211
    for(Component* component : myElements) {
        delete(component);
gsell's avatar
gsell committed
212
    }
213 214
    for(auto fd : FieldDimensions) {
        delete(fd);
gsell's avatar
gsell committed
215 216
    }
    delete itsBeamline;
217
    // delete opalRing_m;
gsell's avatar
gsell committed
218 219
}

220 221 222 223 224 225
/**
 * AAA
 *
 * @param none
 */
void ParallelCyclotronTracker::initializeBoundaryGeometry() {
226 227
    for(Component * component : myElements) {
        bgf_m = dynamic_cast<ElementBase *>(component)->getBoundaryGeometry();
228 229 230 231 232 233 234 235 236 237
        if(!bgf_m)
            continue;
        else
            break;
    }
    if (bgf_m) {
        itsDataSink->writeGeomToVtk(*bgf_m, std::string("data/testGeometry-00000.vtk"));
        OpalData::getInstance()->setGlobalGeometry(bgf_m);
        *gmsg << "* Boundary geometry initialized " << endl;
    }
238
}
239

240 241 242 243 244 245
/**
 *
 *
 * @param fn Base file name
 */
void ParallelCyclotronTracker::bgf_main_collision_test() {
246
    if(!bgf_m) return;
247

248
    Inform msg("bgf_main_collision_test ");
249

250 251 252
    /**
     *Here we check if a particles is outside the domain, flag it for deletion
     */
253

254
    Vector_t intecoords = 0.0;
255

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

259
    int triId = 0;
frey_m's avatar
frey_m committed
260 261 262
    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);
263
        if(res >= 0) {
frey_m's avatar
frey_m committed
264
            itsBunch_m->Bin[i] = -1;
265
        }
266
    }
267 268 269
}


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

277
    std::string  SfileName2 = SfileName + std::string("-Angle0.dat");
gsell's avatar
gsell committed
278
    outfTheta0_m.precision(8);
279
    outfTheta0_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
280
    outfTheta0_m.open(SfileName2.c_str());
281 282 283
    outfTheta0_m << "#  r [mm]      beta_r*gamma       "
                 << "theta [mm]      beta_theta*gamma        "
                 << "z [mm]          beta_z*gamma" << std::endl;
gsell's avatar
gsell committed
284

285
    SfileName2 = SfileName + std::string("-Angle1.dat");
gsell's avatar
gsell committed
286
    outfTheta1_m.precision(8);
287
    outfTheta1_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
288
    outfTheta1_m.open(SfileName2.c_str());
289 290 291
    outfTheta1_m << "#  r [mm]      beta_r*gamma       "
                 << "theta [mm]      beta_theta*gamma        "
                 << "z [mm]          beta_z*gamma"  << std::endl;
gsell's avatar
gsell committed
292

293
    SfileName2 = SfileName + std::string("-Angle2.dat");
gsell's avatar
gsell committed
294
    outfTheta2_m.precision(8);
295
    outfTheta2_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
296
    outfTheta2_m.open(SfileName2.c_str());
297 298 299
    outfTheta2_m << "#  r [mm]      beta_r*gamma       "
                 << "theta [mm]      beta_theta*gamma        "
                 << "z [mm]          beta_z*gamma"  << std::endl;
gsell's avatar
gsell committed
300 301 302

    // for single Particle Mode, output after each turn, to define matched initial phase ellipse.

303
    SfileName2 = SfileName + std::string("-afterEachTurn.dat");
gsell's avatar
gsell committed
304 305

    outfThetaEachTurn_m.precision(8);
306
    outfThetaEachTurn_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
307 308

    outfThetaEachTurn_m.open(SfileName2.c_str());
309 310 311
    outfTheta2_m << "#  r [mm]      beta_r*gamma       "
                 << "theta [mm]      beta_theta*gamma        "
                 << "z [mm]          beta_z*gamma"  << std::endl;
gsell's avatar
gsell committed
312 313 314 315 316 317 318 319 320 321 322 323 324 325 326
}

/**
 * Close all files related to
 * special output in the Cyclotron
 * mode.
 */
void ParallelCyclotronTracker::closeFiles() {

    outfTheta0_m.close();
    outfTheta1_m.close();
    outfTheta2_m.close();
    outfThetaEachTurn_m.close();
}

327
/**
328 329 330
 *
 * @param ring
 */
kraus's avatar
kraus committed
331
void ParallelCyclotronTracker::visitRing(const Ring &ring) {
Daniel Winklehner's avatar
Daniel Winklehner committed
332

333
    *gmsg << "* ----------------------------- Adding Ring ------------------------------ *" << endl;
Daniel Winklehner's avatar
Daniel Winklehner committed
334 335

    if (opalRing_m != NULL)
336
        delete opalRing_m;
Daniel Winklehner's avatar
Daniel Winklehner committed
337

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

340
    myElements.push_back(opalRing_m);
Daniel Winklehner's avatar
Daniel Winklehner committed
341

frey_m's avatar
frey_m committed
342
    opalRing_m->initialise(itsBunch_m);
343 344 345 346

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

348
    if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
kraus's avatar
kraus committed
349
        throw OpalException("Error in ParallelCyclotronTracker::visitRing",
350 351
                            "PHIINIT is out of [-180, 180)!");
    }
Daniel Winklehner's avatar
Daniel Winklehner committed
352 353

    referenceZ = 0.0;
354
    referencePz = 0.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
355 356 357 358

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

359 360
    if(referencePtot < 0.0)
        referencePt *= -1.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
361

362 363
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
gsell's avatar
gsell committed
364

365
    double BcParameter[8];
Daniel Winklehner's avatar
Daniel Winklehner committed
366

367
    for(int i = 0; i < 8; i++)
Daniel Winklehner's avatar
Daniel Winklehner committed
368 369
        BcParameter[i] = 0.0;

370
    buildupFieldList(BcParameter, ElementBase::RING, opalRing_m);
371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386

    // 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
387 388 389 390 391 392 393 394

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

395
    *gmsg << "* -------------------------- Adding Cyclotron ---------------------------- *" << endl;
gsell's avatar
gsell committed
396

397 398
    Cyclotron *elptr = dynamic_cast<Cyclotron *>(cycl.clone());
    myElements.push_back(elptr);
399

400
    // Is this a Spiral Inflector Simulation? If yes, we'll give the user some
401 402 403 404 405 406
    // useful information
    spiral_flag = elptr->getSpiralFlag();

    if(spiral_flag) {

        *gmsg << endl << "* This is a Spiral Inflector Simulation! This means the following:" << endl;
407
        *gmsg         << "* 1.) It is up to the user to provide appropriate geometry, electric and magnetic fields!" << endl;
408 409
        *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;
410
        *gmsg         << "* 2.) For high currentst is strongly recommended to use the SAAMG fieldsolver," << endl;
411 412
        *gmsg         << "*     FFT does not give the correct results (boundaty conditions are missing)." << endl;
        *gmsg         << "* 3.) The whole geometry will be meshed and used for the fieldsolve." << endl;
413 414
        *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;
415 416 417 418 419 420 421
        *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;
        numBunch_m = 1;

    }

422
    // Fresh run (no restart):
423
    if(!OpalData::getInstance()->inRestartRun()) {
424

425 426
        // Get reference values from cyclotron element 
        // For now, these are still stored in mm. should be the only ones. -DW
427 428 429 430
        referenceR     = elptr->getRinit();
        referenceTheta = elptr->getPHIinit();
        referenceZ     = elptr->getZinit();
        referencePr    = elptr->getPRinit();
431 432
        referencePz    = elptr->getPZinit();

433
        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
434 435
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                "PHIINIT is out of [-180, 180)!");
436 437 438 439 440
        }

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

        // Calculate reference azimuthal (tangential) momentum from total-, z- and radial momentum:
441
        float insqrt = referencePtot * referencePtot - \
442
            referencePr * referencePr - referencePz * referencePz;
443 444 445 446 447 448 449 450 451

        if(insqrt < 0) {

            if(insqrt > -1.0e-10) {

	        referencePt = 0.0;

            } else {

452
	        throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
453
                                    "Pt imaginary!");
454 455 456 457 458 459 460
            }

        } else {

            referencePt = sqrt(insqrt);
        }

461
        if(referencePtot < 0.0)
Daniel Winklehner's avatar
Daniel Winklehner committed
462
            referencePt *= -1.0;
463 464
        // End calculate referencePt

465
        // Restart a run:
466 467
    } else {

468
        // If the user wants to save the restarted run in local frame,
469
        // make sure the previous h5 file was local too
470
      if (Options::psDumpFrame != Options::GLOBAL) {
winklehner_d's avatar
-DW  
winklehner_d committed
471
        if (!previousH5Local) {
472 473
                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
474
        }
475 476
            // 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
477 478
      } else {
        if (previousH5Local) {
479 480
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                    "You are trying a global restart from a local h5 file!");
481
            }
482
        }
483

484
        // Adjust some of the reference variables from the h5 file
485 486
        referencePhi *= Physics::deg2rad;
        referencePsi *= Physics::deg2rad;
487
        referencePtot = bega;
488
        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
489
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
490 491
                                "PHIINIT is out of [-180, 180)!");
        }
492 493
    }

494 495
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
496

497
    *gmsg << endl;
adelmann's avatar
adelmann committed
498
    *gmsg << "* Bunch global starting position:" << endl;
499 500
    *gmsg << "* RINIT = " << referenceR  << " [mm]" << endl;
    *gmsg << "* PHIINIT = " << referenceTheta << " [deg]" << endl;
501
    *gmsg << "* ZINIT = " << referenceZ << " [mm]" << endl;
502
    *gmsg << endl;
adelmann's avatar
adelmann committed
503
    *gmsg << "* Bunch global starting momenta:" << endl;
gsell's avatar
gsell committed
504 505
    *gmsg << "* Initial gamma = " << itsReference.getGamma() << endl;
    *gmsg << "* Initial beta = " << itsReference.getBeta() << endl;
506
    *gmsg << "* Reference total momentum (beta * gamma) = " << referencePtot * 1000.0 << " [MCU]" << endl;
507 508 509
    *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;
510
    *gmsg << endl;
adelmann's avatar
adelmann committed
511

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

515 516 517
    // 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
518

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

522
    std::string type = elptr->getCyclotronType();
523
    *gmsg << "* Type of cyclotron = " << type << " " << endl;
524

525 526
    double rmin = elptr->getMinR();
    double rmax = elptr->getMaxR();
adelmann's avatar
adelmann committed
527
    *gmsg << "* Radial aperture = " << rmin << " ... " << rmax<<" [m] "<< endl;
528 529 530

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

    double h = elptr->getCyclHarm();
534
    *gmsg << "* Number of trimcoils = " << elptr->getNumberOfTrimcoils() << endl;
535
    *gmsg << "* Harmonic number h = " << h << " " << endl;
536

gsell's avatar
gsell committed
537 538 539 540 541 542 543 544 545
    /**
     * 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
     * fieldflag = 5, readin FFAG format file for MSU/FNAL FFAG
     * fieldflag = 6, readin both median plane B field map and 3D E field map of RF cavity for compact cyclotron
546
     * fieldflag = 7, read in fields for Daniel's synchrocyclotron simulations
gsell's avatar
gsell committed
547 548
     */
    int  fieldflag;
549
    if(type == std::string("CARBONCYCL")) {
gsell's avatar
gsell committed
550
        fieldflag = 2;
551
    } else if(type == std::string("CYCIAE")) {
gsell's avatar
gsell committed
552
        fieldflag = 3;
553
    } else if(type == std::string("AVFEQ")) {
gsell's avatar
gsell committed
554
        fieldflag = 4;
555
    } else if(type == std::string("FFAG")) {
gsell's avatar
gsell committed
556
        fieldflag = 5;
557
    } else if(type == std::string("BANDRF")) {
gsell's avatar
gsell committed
558
        fieldflag = 6;
559 560
    } else if(type == std::string("SYNCHROCYCLOTRON")) {
	fieldflag = 7;
561
    } else //(type == "RING")
gsell's avatar
gsell committed
562 563
        fieldflag = 1;

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

    double BcParameter[8];
568

569
    for(int i = 0; i < 8; i++)
570
        BcParameter[i] = 0.0;
571

572 573
    BcParameter[0] = 0.001 * elptr->getRmin();
    BcParameter[1] = 0.001 * elptr->getRmax();
gsell's avatar
gsell committed
574

575
    // Store inner radius and outer radius of cyclotron field map in the list
576
    buildupFieldList(BcParameter, ElementBase::CYCLOTRON, elptr);
gsell's avatar
gsell committed
577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593
}

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

/**
 *
 *
 * @param coll
 */
void ParallelCyclotronTracker::visitCollimator(const Collimator &coll) {

594
    *gmsg << "* --------- Collimator -----------------------------" << endl;
gsell's avatar
gsell committed
595

596 597
    Collimator* elptr = dynamic_cast<Collimator *>(coll.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
598

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

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

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

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

611
    double zstart = elptr->getZStart();
adelmann's avatar
adelmann committed
612
    *gmsg << "* Zstart= " << zstart << " [mm]" << endl;
613 614

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

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

frey_m's avatar
frey_m committed
620
    elptr->initialise(itsBunch_m);
gsell's avatar
gsell committed
621 622 623 624

    double BcParameter[8];
    for(int i = 0; i < 8; i++)
        BcParameter[i] = 0.0;
625

626 627 628 629 630 631
    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 ;

632
    buildupFieldList(BcParameter, ElementBase::COLLIMATOR, elptr);
gsell's avatar
gsell committed
633 634 635 636 637 638 639 640 641 642 643 644
}

/**
 *
 *
 * @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
645 646 647 648 649 650 651 652 653 654 655 656
/**
 *
 *
 * @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
657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686
/**
 *
 *
 * @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()));
}

/**
 *
 *
 * @param lamb
 */
void ParallelCyclotronTracker::visitLambertson(const Lambertson &lamb) {
    *gmsg << "In Lambertson; L= " << lamb.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Lambertson *>(lamb.clone()));
}

687 688 689
void ParallelCyclotronTracker::visitOffset(const Offset & off) {
    if (opalRing_m == NULL)
        throw OpalException(
690 691
                            "ParallelCylcotronTracker::visitOffset",
                            "Attempt to place an offset when Ring not defined");
692 693
    Offset* offNonConst = const_cast<Offset*>(&off);
    offNonConst->updateGeometry(opalRing_m->getNextPosition(),
694
                                opalRing_m->getNextNormal());
695 696 697
    opalRing_m->appendElement(off);
}

gsell's avatar
gsell committed
698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719
/**
 *
 *
 * @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());
}

720

gsell's avatar
gsell committed
721 722 723 724 725 726 727 728 729 730
/**
 *
 *
 * @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
731 732 733 734 735 736 737
/**
 *
 *
 * @param multT
 */
void ParallelCyclotronTracker::visitMultipoleT(const MultipoleT &multT) {
    *gmsg << "Adding MultipoleT" << endl;
738
    if (opalRing_m != NULL) {
ext-rogers_c's avatar
ext-rogers_c committed
739
        opalRing_m->appendElement(multT);
740
    } else {
ext-rogers_c's avatar
ext-rogers_c committed
741 742
        throw OpalException("ParallelCyclotronTracker::visitMultipoleT",
                            "Need to define a RINGDEFINITION to use MultipoleT element");
743
    }
ext-rogers_c's avatar
ext-rogers_c committed
744 745 746
    myElements.push_back(dynamic_cast<MultipoleT *>(multT.clone()));
}

gsell's avatar
gsell committed
747 748 749 750 751 752
/**
 *
 *
 * @param prob
 */
void ParallelCyclotronTracker::visitProbe(const Probe &prob) {
753
    *gmsg << "* -----------  Probe -------------------------------" << endl;
754 755
    Probe *elptr = dynamic_cast<Probe *>(prob.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
756

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

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

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

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

769
    double width = elptr->getWidth();
770
    *gmsg << "Width= " << width << " [mm]" << endl;
gsell's avatar
gsell committed
771 772 773


    // initialise, do nothing
frey_m's avatar
frey_m committed
774
    elptr->initialise(itsBunch_m);
gsell's avatar
gsell committed
775 776 777 778

    double BcParameter[8];
    for(int i = 0; i < 8; i++)
        BcParameter[i] = 0.0;
779

780 781 782 783 784
    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 ;
gsell's avatar
gsell committed
785 786

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

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

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

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

819 820 821 822 823 824
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",
825
                            "Need to define a RINGDEFINITION to use VariableRFCavity element");
826 827
}

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

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

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

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

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

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

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

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

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

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

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

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

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

875 876 877 878
    std::shared_ptr<AbstractTimeDependence> freq_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> ampl_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> phase_atd = nullptr;

879
    dvector_t  unityVec;
880 881 882 883
    unityVec.push_back(1.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
884

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

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

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

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