ParallelCyclotronTracker.cpp 124 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
#include "AbsBeamline/CCollimator.h"
gsell's avatar
gsell committed
31 32
#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
#include "AbsBeamline/VariableRFCavityFringeField.h"
57

58 59
#include "AbstractObjects/Element.h"

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

#include "BeamlineGeometry/Euclid3D.h"
#include "BeamlineGeometry/PlanarArcGeometry.h"
Jianjun Yang's avatar
Jianjun Yang committed
66
#include "BeamlineGeometry/RBendGeometry.h"
ext-rogers_c's avatar
ext-rogers_c committed
67
#include "BeamlineGeometry/StraightGeometry.h"
gsell's avatar
gsell committed
68 69 70 71 72 73 74 75 76 77 78
#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/OpalException.h"

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

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

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

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

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

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

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

104

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

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

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

extern Inform *gmsg;

// typedef FVector<double, PSdim> Vector;

/**
 * Constructor ParallelCyclotronTracker
 *
 * @param beamline
 * @param reference
 * @param revBeam
 * @param revTrack
 */
ParallelCyclotronTracker::ParallelCyclotronTracker(const Beamline &beamline,
126 127
                                                   const PartData &reference,
                                                   bool revBeam, bool revTrack):
gsell's avatar
gsell committed
128
    Tracker(beamline, reference, revBeam, revTrack),
129 130
    itsDataSink(nullptr),
    bgf_m(nullptr),
131
    lastDumpedStep_m(0),
132
    eta_m(0.01),
gsell's avatar
gsell committed
133 134
    myNode_m(Ippl::myNode()),
    initialLocalNum_m(0),
135
    initialTotalNum_m(0),
136
    onebunch_m(OpalData::getInstance()->getInputBasename() + "-onebunch.h5"),
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),
gsell's avatar
gsell committed
167
    myNode_m(Ippl::myNode()),
frey_m's avatar
frey_m committed
168 169
    initialLocalNum_m(bunch->getLocalNum()),
    initialTotalNum_m(bunch->getTotalNum()),
170
    onebunch_m(OpalData::getInstance()->getInputBasename() + "-onebunch.h5"),
171 172
    opalRing_m(nullptr),
    itsStepper_mp(nullptr) {
gsell's avatar
gsell committed
173 174
    itsBeamline = dynamic_cast<Beamline *>(beamline.clone());
    itsDataSink = &ds;
frey_m's avatar
frey_m committed
175
    //  scaleFactor_m = itsBunch_m->getdT() * c;
gsell's avatar
gsell committed
176
    scaleFactor_m = 1;
177
    multiBunchMode_m = MB_MODE::NONE;
gsell's avatar
gsell committed
178 179 180 181 182

    IntegrationTimer_m = IpplTimings::getTimer("Integration");
    TransformTimer_m   = IpplTimings::getTimer("Frametransform");
    DumpTimer_m        = IpplTimings::getTimer("Dump");
    BinRepartTimer_m   = IpplTimings::getTimer("Binaryrepart");
183

184 185 186 187 188 189 190 191 192
    // 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;
193

194 195 196 197 198 199 200 201
    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
202 203 204 205 206 207 208
}

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

219 220 221 222 223 224
/**
 * AAA
 *
 * @param none
 */
void ParallelCyclotronTracker::initializeBoundaryGeometry() {
225 226
    for(Component * component : myElements) {
        bgf_m = dynamic_cast<ElementBase *>(component)->getBoundaryGeometry();
227 228 229 230 231 232 233 234 235 236
        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;
    }
237
}
238

239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260

/// set the working sub-mode for multi-bunch mode: "FORCE" or "AUTO"
inline
void ParallelCyclotronTracker::setMultiBunchMode(const std::string& mbmode)
{
    if ( mbmode.compare("FORCE") == 0 ) {
        *gmsg << "FORCE mode: The multi bunches will be injected consecutively "
              << "after each revolution, until get \"TURNS\" bunches." << endl;
        multiBunchMode_m = MB_MODE::FORCE;
    } else if ( mbmode.compare("AUTO") == 0 ) {
        *gmsg << "AUTO mode: The multi bunches will be injected only when "
                      << "the distance between two neighboring bunches " << endl
                      << "is below the limitation. The control parameter is set to "
                      << CoeffDBunches_m << endl;
        multiBunchMode_m = MB_MODE::AUTO;
    } else if ( mbmode.compare("NONE") == 0 )
        multiBunchMode_m = MB_MODE::NONE;
    else
        throw OpalException("ParallelCyclotronTracker::setMultiBunchMode()",
                            "MBMODE name \"" + mbmode + "\" unknown.");
}

261 262 263 264 265 266
/**
 *
 *
 * @param fn Base file name
 */
void ParallelCyclotronTracker::bgf_main_collision_test() {
267
    if(!bgf_m) return;
268

269
    Inform msg("bgf_main_collision_test ");
270

271 272 273
    /**
     *Here we check if a particles is outside the domain, flag it for deletion
     */
274

275
    Vector_t intecoords = 0.0;
276

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

280
    int triId = 0;
frey_m's avatar
frey_m committed
281 282 283
    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);
284
        if(res >= 0) {
frey_m's avatar
frey_m committed
285
            itsBunch_m->Bin[i] = -1;
286
        }
287
    }
288 289 290
}


gsell's avatar
gsell committed
291 292 293 294 295
/**
 *
 *
 * @param fn Base file name
 */
296
void ParallelCyclotronTracker::openFiles(std::string SfileName) {
gsell's avatar
gsell committed
297

298
    std::string  SfileName2 = SfileName + std::string("-Angle0.dat");
gsell's avatar
gsell committed
299
    outfTheta0_m.precision(8);
300
    outfTheta0_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
301
    outfTheta0_m.open(SfileName2.c_str());
302 303 304
    outfTheta0_m << "#  r [mm]      beta_r*gamma       "
                 << "theta [mm]      beta_theta*gamma        "
                 << "z [mm]          beta_z*gamma" << std::endl;
gsell's avatar
gsell committed
305

306
    SfileName2 = SfileName + std::string("-Angle1.dat");
gsell's avatar
gsell committed
307
    outfTheta1_m.precision(8);
308
    outfTheta1_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
309
    outfTheta1_m.open(SfileName2.c_str());
310 311 312
    outfTheta1_m << "#  r [mm]      beta_r*gamma       "
                 << "theta [mm]      beta_theta*gamma        "
                 << "z [mm]          beta_z*gamma"  << std::endl;
gsell's avatar
gsell committed
313

314
    SfileName2 = SfileName + std::string("-Angle2.dat");
gsell's avatar
gsell committed
315
    outfTheta2_m.precision(8);
316
    outfTheta2_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
317
    outfTheta2_m.open(SfileName2.c_str());
318 319 320
    outfTheta2_m << "#  r [mm]      beta_r*gamma       "
                 << "theta [mm]      beta_theta*gamma        "
                 << "z [mm]          beta_z*gamma"  << std::endl;
gsell's avatar
gsell committed
321 322 323

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

324
    SfileName2 = SfileName + std::string("-afterEachTurn.dat");
gsell's avatar
gsell committed
325 326

    outfThetaEachTurn_m.precision(8);
327
    outfThetaEachTurn_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
328 329

    outfThetaEachTurn_m.open(SfileName2.c_str());
330 331 332
    outfTheta2_m << "#  r [mm]      beta_r*gamma       "
                 << "theta [mm]      beta_theta*gamma        "
                 << "z [mm]          beta_z*gamma"  << std::endl;
gsell's avatar
gsell committed
333 334 335 336 337 338 339 340 341 342 343 344 345 346 347
}

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

348
/**
349 350 351
 *
 * @param ring
 */
kraus's avatar
kraus committed
352
void ParallelCyclotronTracker::visitRing(const Ring &ring) {
Daniel Winklehner's avatar
Daniel Winklehner committed
353

354
    *gmsg << "* ----------------------------- Adding Ring ------------------------------ *" << endl;
Daniel Winklehner's avatar
Daniel Winklehner committed
355

snuverink_j's avatar
snuverink_j committed
356
    delete opalRing_m;
Daniel Winklehner's avatar
Daniel Winklehner committed
357

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

360
    myElements.push_back(opalRing_m);
Daniel Winklehner's avatar
Daniel Winklehner committed
361

frey_m's avatar
frey_m committed
362
    opalRing_m->initialise(itsBunch_m);
363 364 365 366

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

368
    if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
kraus's avatar
kraus committed
369
        throw OpalException("Error in ParallelCyclotronTracker::visitRing",
370 371
                            "PHIINIT is out of [-180, 180)!");
    }
Daniel Winklehner's avatar
Daniel Winklehner committed
372 373

    referenceZ = 0.0;
374
    referencePz = 0.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
375 376 377 378

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

379 380
    if(referencePtot < 0.0)
        referencePt *= -1.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
381

382 383
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
gsell's avatar
gsell committed
384

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

387
    buildupFieldList(BcParameter, ElementBase::RING, opalRing_m);
388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403

    // 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
404 405 406 407 408 409 410 411

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

412
    *gmsg << "* -------------------------- Adding Cyclotron ---------------------------- *" << endl;
gsell's avatar
gsell committed
413

414 415
    Cyclotron *elptr = dynamic_cast<Cyclotron *>(cycl.clone());
    myElements.push_back(elptr);
416

417
    // Is this a Spiral Inflector Simulation? If yes, we'll give the user some
418 419 420 421 422 423
    // useful information
    spiral_flag = elptr->getSpiralFlag();

    if(spiral_flag) {

        *gmsg << endl << "* This is a Spiral Inflector Simulation! This means the following:" << endl;
424
        *gmsg         << "* 1.) It is up to the user to provide appropriate geometry, electric and magnetic fields!" << endl;
425 426
        *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
427
        *gmsg         << "* 2.) For high currents it is strongly recommended to use the SAAMG fieldsolver," << endl;
428 429
        *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;
430 431
        *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;
432 433 434 435 436 437 438
        *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;

    }

439
    // Fresh run (no restart):
440
    if(!OpalData::getInstance()->inRestartRun()) {
441

442
        // Get reference values from cyclotron element
443
        // For now, these are still stored in mm. should be the only ones. -DW
444 445 446 447
        referenceR     = elptr->getRinit();
        referenceTheta = elptr->getPHIinit();
        referenceZ     = elptr->getZinit();
        referencePr    = elptr->getPRinit();
448 449
        referencePz    = elptr->getPZinit();

450
        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
451 452
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                "PHIINIT is out of [-180, 180)!");
453 454 455 456 457
        }

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

        // Calculate reference azimuthal (tangential) momentum from total-, z- and radial momentum:
458
        float insqrt = referencePtot * referencePtot - \
459
            referencePr * referencePr - referencePz * referencePz;
460 461 462 463 464 465 466 467 468

        if(insqrt < 0) {

            if(insqrt > -1.0e-10) {

	        referencePt = 0.0;

            } else {

469
	        throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
470
                                    "Pt imaginary!");
471 472 473 474 475 476 477
            }

        } else {

            referencePt = sqrt(insqrt);
        }

478
        if(referencePtot < 0.0)
Daniel Winklehner's avatar
Daniel Winklehner committed
479
            referencePt *= -1.0;
480 481
        // End calculate referencePt

482
        // Restart a run:
483 484
    } else {

485
        // If the user wants to save the restarted run in local frame,
486
        // make sure the previous h5 file was local too
487
      if (Options::psDumpFrame != Options::GLOBAL) {
winklehner_d's avatar
-DW  
winklehner_d committed
488
        if (!previousH5Local) {
489 490
                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
491
        }
492 493
            // 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
494 495
      } else {
        if (previousH5Local) {
496 497
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                    "You are trying a global restart from a local h5 file!");
498
            }
499
        }
500

501
        // Adjust some of the reference variables from the h5 file
502 503
        referencePhi *= Physics::deg2rad;
        referencePsi *= Physics::deg2rad;
504
        referencePtot = bega;
505
        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
506
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
507 508
                                "PHIINIT is out of [-180, 180)!");
        }
509 510
    }

511 512
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
513

514
    *gmsg << endl;
adelmann's avatar
adelmann committed
515
    *gmsg << "* Bunch global starting position:" << endl;
516 517
    *gmsg << "* RINIT = " << referenceR  << " [mm]" << endl;
    *gmsg << "* PHIINIT = " << referenceTheta << " [deg]" << endl;
518
    *gmsg << "* ZINIT = " << referenceZ << " [mm]" << endl;
519
    *gmsg << endl;
adelmann's avatar
adelmann committed
520
    *gmsg << "* Bunch global starting momenta:" << endl;
gsell's avatar
gsell committed
521 522
    *gmsg << "* Initial gamma = " << itsReference.getGamma() << endl;
    *gmsg << "* Initial beta = " << itsReference.getBeta() << endl;
523
    *gmsg << "* Reference total momentum (beta * gamma) = " << referencePtot * 1000.0 << " [MCU]" << endl;
524 525 526
    *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;
527
    *gmsg << endl;
adelmann's avatar
adelmann committed
528

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

532 533 534
    // 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
535

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

539
    std::string type = elptr->getCyclotronType();
540
    *gmsg << "* Type of cyclotron = " << type << " " << endl;
541

542 543
    double rmin = elptr->getMinR();
    double rmax = elptr->getMaxR();
adelmann's avatar
adelmann committed
544
    *gmsg << "* Radial aperture = " << rmin << " ... " << rmax<<" [m] "<< endl;
545 546 547

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

    double h = elptr->getCyclHarm();
551
    *gmsg << "* Number of trimcoils = " << elptr->getNumberOfTrimcoils() << endl;
552
    *gmsg << "* Harmonic number h = " << h << " " << endl;
553

gsell's avatar
gsell committed
554 555 556 557 558 559 560 561 562
    /**
     * 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
563
     * fieldflag = 7, read in fields for Daniel's synchrocyclotron simulations
gsell's avatar
gsell committed
564
     */
565
    int fieldflag = elptr->getFieldFlag(type);
gsell's avatar
gsell committed
566

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

snuverink_j's avatar
snuverink_j committed
570
    double BcParameter[8] = {};
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
}

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

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

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

596
    CCollimator* elptr = dynamic_cast<CCollimator *>(coll.clone());
597
    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

snuverink_j's avatar
snuverink_j committed
622
    double BcParameter[8] = {};
623

624 625 626 627 628 629
    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 ;

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

/**
 *
 *
 * @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
643 644 645 646 647 648 649 650 651 652 653 654
/**
 *
 *
 * @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
655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674
/**
 *
 *
 * @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()));
}

675 676
/**
 *
677 678
 *
 *  @param
679 680 681 682 683
 */
void ParallelCyclotronTracker::visitFlexibleCollimator(const FlexibleCollimator &) {

}

gsell's avatar
gsell committed
684 685 686 687 688 689 690 691 692 693
/**
 *
 *
 * @param lamb
 */
void ParallelCyclotronTracker::visitLambertson(const Lambertson &lamb) {
    *gmsg << "In Lambertson; L= " << lamb.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Lambertson *>(lamb.clone()));
}

694 695 696
void ParallelCyclotronTracker::visitOffset(const Offset & off) {
    if (opalRing_m == NULL)
        throw OpalException(
697 698
                            "ParallelCylcotronTracker::visitOffset",
                            "Attempt to place an offset when Ring not defined");
699 700
    Offset* offNonConst = const_cast<Offset*>(&off);
    offNonConst->updateGeometry(opalRing_m->getNextPosition(),
701
                                opalRing_m->getNextNormal());
702 703 704
    opalRing_m->appendElement(off);
}

gsell's avatar
gsell committed
705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726
/**
 *
 *
 * @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());
}

727

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

gsell's avatar
gsell committed
754 755 756 757 758 759
/**
 *
 *
 * @param prob
 */
void ParallelCyclotronTracker::visitProbe(const Probe &prob) {
760
    *gmsg << "* -----------  Probe -------------------------------" << endl;
761 762
    Probe *elptr = dynamic_cast<Probe *>(prob.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
763

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

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

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

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

776
    double width = elptr->getWidth();
777
    *gmsg << "Width= " << width << " [mm]" << endl;
gsell's avatar
gsell committed
778 779 780


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

snuverink_j's avatar
snuverink_j committed
783
    double BcParameter[8] = {};
784

785 786 787 788 789
    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
790 791

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

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

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

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

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

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

gsell's avatar
gsell committed
843 844 845 846 847 848 849 850
/**
 *
 *
 * @param as
 */
void ParallelCyclotronTracker::visitRFCavity(const RFCavity &as) {

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

852 853
    RFCavity *elptr = dynamic_cast<RFCavity *>(as.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
854 855 856 857 858 859 860

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

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

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

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

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

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

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

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

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

885 886 887
    /*
      Setup time dependence and in case of no
      timedependence use a polynom with  a_0 = 1 and a_k = 0, k = 1,2,3.
888
    */
889

890 891 892 893
    std::shared_ptr<AbstractTimeDependence> freq_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> ampl_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> phase_atd = nullptr;

894
    dvector_t  unityVec;
895 896 897 898
    unityVec.push_back(1.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
899

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

    if (elptr->getAmplitudeModelName() != "") {
908 909
        ampl_atd = AbstractTimeDependence::getTimeDependence(elptr->getAmplitudeModelName());
        *gmsg << "* Variable amplitude RF Model name " << elptr->getAmplitudeModelName() << endl;
910
    }
911 912
    else
        ampl_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));