ParallelCyclotronTracker.cpp 126 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 "Beamlines/FlaggedBeamline.h"
60
#include "Elements/OpalBeamline.h"
kraus's avatar
kraus committed
61
#include "AbsBeamline/Ring.h"
gsell's avatar
gsell committed
62 63 64

#include "BeamlineGeometry/Euclid3D.h"
#include "BeamlineGeometry/PlanarArcGeometry.h"
Jianjun Yang's avatar
Jianjun Yang committed
65
#include "BeamlineGeometry/RBendGeometry.h"
ext-rogers_c's avatar
ext-rogers_c committed
66
#include "BeamlineGeometry/StraightGeometry.h"
gsell's avatar
gsell committed
67 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/NumToStr.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),
133 134
    initialR_m(nullptr),
    initialP_m(nullptr),
gsell's avatar
gsell committed
135 136
    myNode_m(Ippl::myNode()),
    initialLocalNum_m(0),
137
    initialTotalNum_m(0),
138
    opalRing_m(nullptr),
139 140
    itsStepper_mp(nullptr),
    mode_m(MODE::UNDEFINED),
141
    stepper_m(stepper::INTEGRATOR::UNDEFINED) {
gsell's avatar
gsell committed
142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157
    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
158
                                                   PartBunchBase<double, 3> *bunch,
gsell's avatar
gsell committed
159 160 161
                                                   DataSink &ds,
                                                   const PartData &reference,
                                                   bool revBeam, bool revTrack,
162
                                                   int maxSTEPS, int timeIntegrator):
163
    Tracker(beamline, bunch, reference, revBeam, revTrack),
164
    bgf_m(nullptr),
gsell's avatar
gsell committed
165
    maxSteps_m(maxSTEPS),
166
    lastDumpedStep_m(0),
167
    eta_m(0.01),
168 169
    initialR_m(nullptr),
    initialP_m(nullptr),
gsell's avatar
gsell committed
170
    myNode_m(Ippl::myNode()),
frey_m's avatar
frey_m committed
171 172
    initialLocalNum_m(bunch->getLocalNum()),
    initialTotalNum_m(bunch->getTotalNum()),
173 174
    opalRing_m(nullptr),
    itsStepper_mp(nullptr) {
gsell's avatar
gsell committed
175 176
    itsBeamline = dynamic_cast<Beamline *>(beamline.clone());
    itsDataSink = &ds;
frey_m's avatar
frey_m committed
177
    //  scaleFactor_m = itsBunch_m->getdT() * c;
gsell's avatar
gsell committed
178 179 180 181 182 183 184
    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");
185

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

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

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

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

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

249
    Inform msg("bgf_main_collision_test ");
250

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

255
    Vector_t intecoords = 0.0;
256

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

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


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

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

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

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

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

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

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

    outfThetaEachTurn_m.open(SfileName2.c_str());
310 311 312
    outfTheta2_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 315 316 317 318 319 320 321 322 323 324 325 326 327
}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    if(spiral_flag) {

        *gmsg << endl << "* This is a Spiral Inflector Simulation! This means the following:" << endl;
408
        *gmsg         << "* 1.) It is up to the user to provide appropriate geometry, electric and magnetic fields!" << endl;
409 410
        *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;
411
        *gmsg         << "* 2.) For high currentst is strongly recommended to use the SAAMG fieldsolver," << endl;
412 413
        *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;
414 415
        *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;
416 417 418 419 420 421 422
        *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;

    }

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

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

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

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

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

        if(insqrt < 0) {

            if(insqrt > -1.0e-10) {

	        referencePt = 0.0;

            } else {

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

        } else {

            referencePt = sqrt(insqrt);
        }

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

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

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

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

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

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

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

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

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

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

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

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

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

gsell's avatar
gsell committed
538 539 540 541 542 543 544 545 546
    /**
     * 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
547
     * fieldflag = 7, read in fields for Daniel's synchrocyclotron simulations
gsell's avatar
gsell committed
548 549
     */
    int  fieldflag;
550
    if(type == std::string("CARBONCYCL")) {
gsell's avatar
gsell committed
551
        fieldflag = 2;
552
    } else if(type == std::string("CYCIAE")) {
gsell's avatar
gsell committed
553
        fieldflag = 3;
554
    } else if(type == std::string("AVFEQ")) {
gsell's avatar
gsell committed
555
        fieldflag = 4;
556
    } else if(type == std::string("FFAG")) {
gsell's avatar
gsell committed
557
        fieldflag = 5;
558
    } else if(type == std::string("BANDRF")) {
gsell's avatar
gsell committed
559
        fieldflag = 6;
560 561
    } else if(type == std::string("SYNCHROCYCLOTRON")) {
	fieldflag = 7;
562
    } else //(type == "RING")
gsell's avatar
gsell committed
563 564
        fieldflag = 1;

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

    double BcParameter[8];
569

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

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

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

/**
 * 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) {

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

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

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

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

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

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

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

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

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

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

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

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

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

/**
 *
 *
 * @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
646 647 648 649 650 651 652 653 654 655 656 657
/**
 *
 *
 * @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
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 687
/**
 *
 *
 * @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()));
}

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

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

721

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

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

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

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

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

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

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


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

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

781 782 783 784 785
    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
786 787

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    if (elptr->getPhaseModelName() != "") {
901 902
        phase_atd = AbstractTimeDependence::getTimeDependence(elptr->getPhaseModelName());
        *gmsg << "* Variable phase RF Model name " << elptr->getPhaseModelName() << endl;
903
    }