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

#include "Algorithms/ParallelCyclotronTracker.h"
20

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

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

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

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

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

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

#include "Physics/Physics.h"

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

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

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

frey_m's avatar
frey_m committed
83 84
//FIXME Remove headers and dynamic_cast in readOneBunchFromFile
#include "Algorithms/PartBunch.h"
85
#ifdef ENABLE_AMR
frey_m's avatar
frey_m committed
86 87
    #include "Algorithms/AmrPartBunch.h"
#endif
frey_m's avatar
frey_m committed
88

gsell's avatar
gsell committed
89 90 91
using Physics::pi;
using Physics::q_e;

92

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

95 96 97
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
98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113

extern Inform *gmsg;

/**
 * Constructor ParallelCyclotronTracker
 *
 * @param beamline
 * @param bunch
 * @param ds
 * @param reference
 * @param revBeam
 * @param revTrack
 * @param maxSTEPS
 * @param timeIntegrator
 */
ParallelCyclotronTracker::ParallelCyclotronTracker(const Beamline &beamline,
frey_m's avatar
frey_m committed
114
                                                   PartBunchBase<double, 3> *bunch,
gsell's avatar
gsell committed
115 116 117
                                                   DataSink &ds,
                                                   const PartData &reference,
                                                   bool revBeam, bool revTrack,
frey_m's avatar
frey_m committed
118
                                                   int maxSTEPS, int timeIntegrator,
frey_m's avatar
frey_m committed
119
                                                   int numBunch)
frey_m's avatar
frey_m committed
120 121 122
    : Tracker(beamline, bunch, reference, revBeam, revTrack)
    , bgf_m(nullptr)
    , maxSteps_m(maxSTEPS)
123
    , numBunch_m(numBunch)
frey_m's avatar
frey_m committed
124 125 126 127 128 129 130 131 132
    , lastDumpedStep_m(0)
    , eta_m(0.01)
    , myNode_m(Ippl::myNode())
    , initialLocalNum_m(bunch->getLocalNum())
    , initialTotalNum_m(bunch->getTotalNum())
    , onebunch_m(OpalData::getInstance()->getInputBasename() + "-onebunch.h5")
    , opalRing_m(nullptr)
    , itsStepper_mp(nullptr)
{
gsell's avatar
gsell committed
133 134
    itsBeamline = dynamic_cast<Beamline *>(beamline.clone());
    itsDataSink = &ds;
135
    multiBunchMode_m = MB_MODE::NONE;
gsell's avatar
gsell committed
136 137 138 139 140

    IntegrationTimer_m = IpplTimings::getTimer("Integration");
    TransformTimer_m   = IpplTimings::getTimer("Frametransform");
    DumpTimer_m        = IpplTimings::getTimer("Dump");
    BinRepartTimer_m   = IpplTimings::getTimer("Binaryrepart");
141 142
    PluginElemTimer_m  = IpplTimings::getTimer("PluginElements");
    DelParticleTimer_m = IpplTimings::getTimer("DeleteParticles");
143

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

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

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

179 180 181 182
/// set the working sub-mode for multi-bunch mode: "FORCE" or "AUTO"
void ParallelCyclotronTracker::setMultiBunchMode(const std::string& mbmode)
{
    if ( mbmode.compare("FORCE") == 0 ) {
183 184
        *gmsg << "FORCE mode: The multi bunches will be injected consecutively" << endl
              << "            after each revolution, until get \"TURNS\" bunches." << endl;
185 186
        multiBunchMode_m = MB_MODE::FORCE;
    } else if ( mbmode.compare("AUTO") == 0 ) {
187 188 189 190
        *gmsg << "AUTO mode: The multi bunches will be injected only when the" << endl
              << "           distance between two neighboring bunches is below" << endl
              << "           the limitation. The control parameter is set to "
              << CoeffDBunches_m << endl;
191 192 193 194 195 196 197 198
        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.");
}

frey_m's avatar
frey_m committed
199
void ParallelCyclotronTracker::setMultiBunchBinning(std::string binning) {
200

frey_m's avatar
frey_m committed
201
    binning = Util::toUpper(binning);
202

frey_m's avatar
frey_m committed
203 204 205 206 207 208 209 210 211 212 213 214
    if ( binning.compare("BUNCH") == 0 ) {
        *gmsg << "Use 'BUNCH' injection for binnning." << endl;
        binningType_m = MB_BINNING::BUNCH;
    } else if ( binning.compare("GAMMA") == 0 ) {
        *gmsg << "Use 'GAMMA' for binning." << endl;
        binningType_m = MB_BINNING::GAMMA;
    } else {
        throw OpalException("ParallelCyclotronTracker::setMultiBunchBinning()",
                            "MB_BINNING name \"" + binning + "\" unknown.");
    }
}

215 216 217 218 219 220
/**
 *
 *
 * @param fn Base file name
 */
void ParallelCyclotronTracker::bgf_main_collision_test() {
221
    if(!bgf_m) return;
222

223
    Inform msg("bgf_main_collision_test ");
224

225
    /**
226
     *Here we check if a particle is outside the domain, flag it for deletion
227
     */
228

229
    Vector_t intecoords = 0.0;
230

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

234
    int triId = 0;
frey_m's avatar
frey_m committed
235 236 237
    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);
238
        if(res >= 0) {
frey_m's avatar
frey_m committed
239
            itsBunch_m->Bin[i] = -1;
240
        }
241
    }
242 243
}

244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261
// only used for dumping into stat file
void ParallelCyclotronTracker::dumpAngle_m(const double& theta) {
    if ( prevAzimuth_m < 0.0 ) { // only at first occurrence
        azimuth_m = theta;
    } else {
        double dtheta = theta - prevAzimuth_m;
        if ( dtheta < 0 ) {
            dtheta += 360.0;
        }
        if ( dtheta > 180 ) { // rotating clockwise, reduce angle
            dtheta -= 360;
        }
        azimuth_m += dtheta;
    }
    prevAzimuth_m = theta;
}


262

gsell's avatar
gsell committed
263 264 265 266 267
/**
 *
 *
 * @param fn Base file name
 */
268
void ParallelCyclotronTracker::openFiles(size_t numFiles, std::string SfileName) {
gsell's avatar
gsell committed
269

270
    for (unsigned int i=0; i<numFiles; i++) {
271 272 273 274 275 276 277 278
        std::string SfileName2 = SfileName;
        if (i<=2) {
            SfileName2 += std::string("-Angle" + std::to_string(int(i)) + ".dat");
        }
        else {
            // for single Particle Mode, output after each turn, to define matched initial phase ellipse.
            SfileName2 += std::string("-afterEachTurn.dat");
        }
gsell's avatar
gsell committed
279

280 281 282 283 284 285
        outfTheta_m.emplace_back(new std::ofstream(SfileName2.c_str()));
        outfTheta_m.back()->precision(8);
        outfTheta_m.back()->setf(std::ios::scientific, std::ios::floatfield);
        *outfTheta_m.back() << "# r [mm]        beta_r*gamma       "
                            << "theta [deg]     beta_theta*gamma        "
                            << "z [mm]          beta_z*gamma" << std::endl;
286
    }
gsell's avatar
gsell committed
287 288 289 290 291 292 293 294
}

/**
 * Close all files related to
 * special output in the Cyclotron
 * mode.
 */
void ParallelCyclotronTracker::closeFiles() {
295
    for (auto & file : outfTheta_m) {
296
        file->close();
297
    }
gsell's avatar
gsell committed
298 299
}

300
/**
301 302 303
 *
 * @param ring
 */
kraus's avatar
kraus committed
304
void ParallelCyclotronTracker::visitRing(const Ring &ring) {
Daniel Winklehner's avatar
Daniel Winklehner committed
305

306
    *gmsg << "* ----------------------------- Adding Ring ------------------------------ *" << endl;
Daniel Winklehner's avatar
Daniel Winklehner committed
307

snuverink_j's avatar
snuverink_j committed
308
    delete opalRing_m;
Daniel Winklehner's avatar
Daniel Winklehner committed
309

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

312
    myElements.push_back(opalRing_m);
Daniel Winklehner's avatar
Daniel Winklehner committed
313

frey_m's avatar
frey_m committed
314
    opalRing_m->initialise(itsBunch_m);
315 316 317 318

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

320
    if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
kraus's avatar
kraus committed
321
        throw OpalException("Error in ParallelCyclotronTracker::visitRing",
322 323
                            "PHIINIT is out of [-180, 180)!");
    }
Daniel Winklehner's avatar
Daniel Winklehner committed
324 325

    referenceZ = 0.0;
326
    referencePz = 0.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
327 328 329 330

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

331 332
    if(referencePtot < 0.0)
        referencePt *= -1.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
333

334 335
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
gsell's avatar
gsell committed
336

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

339
    buildupFieldList(BcParameter, ElementBase::RING, opalRing_m);
340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355

    // 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
356 357 358 359 360 361 362 363

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

364
    *gmsg << "* -------------------------- Adding Cyclotron ---------------------------- *" << endl;
gsell's avatar
gsell committed
365

366 367
    Cyclotron *elptr = dynamic_cast<Cyclotron *>(cycl.clone());
    myElements.push_back(elptr);
368

369
    // Is this a Spiral Inflector Simulation? If yes, we'll give the user some
370 371 372 373 374 375
    // useful information
    spiral_flag = elptr->getSpiralFlag();

    if(spiral_flag) {

        *gmsg << endl << "* This is a Spiral Inflector Simulation! This means the following:" << endl;
376
        *gmsg         << "* 1.) It is up to the user to provide appropriate geometry, electric and magnetic fields!" << endl;
377 378
        *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
379
        *gmsg         << "* 2.) For high currents it is strongly recommended to use the SAAMG fieldsolver," << endl;
380 381
        *gmsg         << "*     FFT does not give the correct results (boundary conditions are missing)." << endl;
        *gmsg         << "* 3.) The whole geometry will be meshed and used for the fieldsolver." << endl;
382 383
        *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;
384 385 386 387 388 389 390
        *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;

    }

391
    // Fresh run (no restart):
392
    if(!OpalData::getInstance()->inRestartRun()) {
393

394
        // Get reference values from cyclotron element
395
        // For now, these are still stored in mm. should be the only ones. -DW
396 397 398 399
        referenceR     = elptr->getRinit();
        referenceTheta = elptr->getPHIinit();
        referenceZ     = elptr->getZinit();
        referencePr    = elptr->getPRinit();
400 401
        referencePz    = elptr->getPZinit();

402
        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
403 404
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                "PHIINIT is out of [-180, 180)!");
405 406 407 408 409
        }

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

        // Calculate reference azimuthal (tangential) momentum from total-, z- and radial momentum:
410
        float insqrt = referencePtot * referencePtot - \
411
            referencePr * referencePr - referencePz * referencePz;
412 413 414 415 416

        if(insqrt < 0) {

            if(insqrt > -1.0e-10) {

417
                referencePt = 0.0;
418 419 420

            } else {

421
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
422
                                    "Pt imaginary!");
423 424 425 426 427 428 429
            }

        } else {

            referencePt = sqrt(insqrt);
        }

430
        if(referencePtot < 0.0)
Daniel Winklehner's avatar
Daniel Winklehner committed
431
            referencePt *= -1.0;
432 433
        // End calculate referencePt

434
        // Restart a run:
435 436
    } else {

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

453
        // Adjust some of the reference variables from the h5 file
454 455
        referencePhi *= Physics::deg2rad;
        referencePsi *= Physics::deg2rad;
456
        referencePtot = bega;
457
        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
458
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
459 460
                                "PHIINIT is out of [-180, 180)!");
        }
461 462
    }

463 464
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
465

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

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

484 485 486
    // 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
487

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

491
    std::string type = elptr->getCyclotronType();
492
    *gmsg << "* Type of cyclotron = " << type << " " << endl;
493

494 495
    double rmin = elptr->getMinR();
    double rmax = elptr->getMaxR();
adelmann's avatar
adelmann committed
496
    *gmsg << "* Radial aperture = " << rmin << " ... " << rmax<<" [m] "<< endl;
497 498 499

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

    double h = elptr->getCyclHarm();
503
    *gmsg << "* Number of trimcoils = " << elptr->getNumberOfTrimcoils() << endl;
504
    *gmsg << "* Harmonic number h = " << h << " " << endl;
505

gsell's avatar
gsell committed
506 507 508 509 510 511 512
    /**
     * To ease the initialise() function, set a integral parameter fieldflag internally.
     * Its value is  by the option "TYPE" of the element  "CYCLOTRON"
     * fieldflag = 1, readin PSI format measured field file (default)
     * fieldflag = 2, readin carbon cyclotron field file created by Jianjun Yang, TYPE=CARBONCYCL
     * fieldflag = 3, readin ANSYS format file for CYCIAE-100 created by Jianjun Yang, TYPE=CYCIAE
     * fieldflag = 4, readin AVFEQ format file for Riken cyclotrons
ext-rogers_c's avatar
ext-rogers_c committed
513
     * fieldflag = 5, readin FFA format file for MSU/FNAL FFA
gsell's avatar
gsell committed
514
     * fieldflag = 6, readin both median plane B field map and 3D E field map of RF cavity for compact cyclotron
515
     * fieldflag = 7, read in fields for Daniel's synchrocyclotron simulations
gsell's avatar
gsell committed
516
     */
517
    int fieldflag = elptr->getFieldFlag(type);
gsell's avatar
gsell committed
518

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

snuverink_j's avatar
snuverink_j committed
522
    double BcParameter[8] = {};
523

524 525
    BcParameter[0] = 0.001 * elptr->getRmin();
    BcParameter[1] = 0.001 * elptr->getRmax();
gsell's avatar
gsell committed
526

527
    // Store inner radius and outer radius of cyclotron field map in the list
528
    buildupFieldList(BcParameter, ElementBase::CYCLOTRON, elptr);
gsell's avatar
gsell committed
529 530 531 532 533 534 535 536 537 538 539 540 541 542 543
}

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

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

546
    *gmsg << "* --------- Collimator -----------------------------" << endl;
gsell's avatar
gsell committed
547

548
    CCollimator* elptr = dynamic_cast<CCollimator *>(coll.clone());
549
    myElements.push_back(elptr);
gsell's avatar
gsell committed
550

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

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

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

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

563
    double zstart = elptr->getZStart();
adelmann's avatar
adelmann committed
564
    *gmsg << "* Zstart= " << zstart << " [mm]" << endl;
565 566

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

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

frey_m's avatar
frey_m committed
572
    elptr->initialise(itsBunch_m);
gsell's avatar
gsell committed
573

snuverink_j's avatar
snuverink_j committed
574
    double BcParameter[8] = {};
575

576 577 578 579 580 581
    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 ;

582
    buildupFieldList(BcParameter, ElementBase::CCOLLIMATOR, elptr);
gsell's avatar
gsell committed
583 584 585 586 587 588 589 590 591 592 593 594
}

/**
 *
 *
 * @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
595 596 597 598 599 600 601 602 603 604 605 606
/**
 *
 *
 * @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
607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626
/**
 *
 *
 * @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()));
}

627 628
/**
 *
629 630
 *
 *  @param
631 632 633 634 635
 */
void ParallelCyclotronTracker::visitFlexibleCollimator(const FlexibleCollimator &) {

}

gsell's avatar
gsell committed
636 637 638 639 640 641 642 643 644 645
/**
 *
 *
 * @param lamb
 */
void ParallelCyclotronTracker::visitLambertson(const Lambertson &lamb) {
    *gmsg << "In Lambertson; L= " << lamb.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Lambertson *>(lamb.clone()));
}

646 647 648
void ParallelCyclotronTracker::visitOffset(const Offset & off) {
    if (opalRing_m == NULL)
        throw OpalException(
649 650
                            "ParallelCylcotronTracker::visitOffset",
                            "Attempt to place an offset when Ring not defined");
651 652
    Offset* offNonConst = const_cast<Offset*>(&off);
    offNonConst->updateGeometry(opalRing_m->getNextPosition(),
653
                                opalRing_m->getNextNormal());
654 655 656
    opalRing_m->appendElement(off);
}

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 687 688
/**
 *
 *
 * @param marker
 */
void ParallelCyclotronTracker::visitMarker(const Marker &marker) {
    //   *gmsg << "In Marker; L= " << marker.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Marker *>(marker.clone()));
    // Do nothing.
}

/**
 *
 *
 * @param corr
 */
void ParallelCyclotronTracker::visitMonitor(const Monitor &corr) {
    //   *gmsg << "In Monitor; L= " << corr.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Monitor *>(corr.clone()));
    //   applyDrift(flip_s * corr.getElementLength());
}

/**
 *
 *
 * @param mult
 */
void ParallelCyclotronTracker::visitMultipole(const Multipole &mult) {
    *gmsg << "In Multipole; L= " << mult.getElementLength() << " however the element is missing " << endl;
    myElements.push_back(dynamic_cast<Multipole *>(mult.clone()));
}

ext-rogers_c's avatar
ext-rogers_c committed
689 690 691 692 693 694 695
/**
 *
 *
 * @param multT
 */
void ParallelCyclotronTracker::visitMultipoleT(const MultipoleT &multT) {
    *gmsg << "Adding MultipoleT" << endl;
696
    if (opalRing_m != NULL) {
ext-rogers_c's avatar
ext-rogers_c committed
697
        opalRing_m->appendElement(multT);
698
    } else {
ext-rogers_c's avatar
ext-rogers_c committed
699 700
        throw OpalException("ParallelCyclotronTracker::visitMultipoleT",
                            "Need to define a RINGDEFINITION to use MultipoleT element");
701
    }
ext-rogers_c's avatar
ext-rogers_c committed
702 703 704
    myElements.push_back(dynamic_cast<MultipoleT *>(multT.clone()));
}

705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752
/**
 *
 *
 * @param multTstraight
 */
void ParallelCyclotronTracker::visitMultipoleTStraight(const MultipoleTStraight &multTstraight) {
    *gmsg << "Adding MultipoleTStraight" << endl;
    if (opalRing_m != NULL) {
        opalRing_m->appendElement(multTstraight);
    } else {
        throw OpalException("ParallelCyclotronTracker::visitMultipoleTStraight",
                            "Need to define a RINGDEFINITION to use MultipoleTStraight element");
    }
    myElements.push_back(dynamic_cast<MultipoleTStraight *>(multTstraight.clone()));
}

/**
 *
 *
 * @param multTccurv
 */
void ParallelCyclotronTracker::visitMultipoleTCurvedConstRadius(const MultipoleTCurvedConstRadius &multTccurv) {
    *gmsg << "Adding MultipoleTCurvedConstRadius" << endl;
    if (opalRing_m != NULL) {
        opalRing_m->appendElement(multTccurv);
    } else {
        throw OpalException("ParallelCyclotronTracker::visitMultipoleTCurvedConstRadius",
                            "Need to define a RINGDEFINITION to use MultipoleTCurvedConstRadius element");
    }
    myElements.push_back(dynamic_cast<MultipoleTCurvedConstRadius *>(multTccurv.clone()));
}

/**
 *
 *
 * @param multTvcurv
 */
void ParallelCyclotronTracker::visitMultipoleTCurvedVarRadius(const MultipoleTCurvedVarRadius &multTvcurv) {
    *gmsg << "Adding MultipoleTCurvedVarRadius" << endl;
    if (opalRing_m != NULL) {
        opalRing_m->appendElement(multTvcurv);
    } else {
        throw OpalException("ParallelCyclotronTracker::visitMultipoleTCurvedVarRadius",
                            "Need to define a RINGDEFINITION to use MultipoleTCurvedVarRadius element");
    }
    myElements.push_back(dynamic_cast<MultipoleTCurvedVarRadius *>(multTvcurv.clone()));
}

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

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

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

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

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

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

snuverink_j's avatar
snuverink_j committed
778
    double BcParameter[8] = {};
779

780 781 782 783
    BcParameter[0] = 0.001 * xstart ;
    BcParameter[1] = 0.001 * xend;
    BcParameter[2] = 0.001 * ystart ;
    BcParameter[3] = 0.001 * yend;
784
    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
}

ext-rogers_c's avatar
ext-rogers_c committed
809 810
void ParallelCyclotronTracker::visitScalingFFAMagnet(const ScalingFFAMagnet &bend) {
    *gmsg << "Adding ScalingFFAMagnet" << 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 {
ext-rogers_c's avatar
ext-rogers_c committed
814 815
        throw OpalException("ParallelCyclotronTracker::visitScalingFFAMagnet",
                            "Need to define a RINGDEFINITION to use ScalingFFAMagnet 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 828 829 830 831 832 833 834 835
}

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

gsell's avatar
gsell committed
838 839 840 841 842 843 844 845
/**
 *
 *
 * @param as
 */
void ParallelCyclotronTracker::visitRFCavity(const RFCavity &as) {

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

847 848
    RFCavity *elptr = dynamic_cast<RFCavity *>(as.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
849 850 851 852 853 854 855

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

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

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

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

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

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

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

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

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

880 881 882
    /*
      Setup time dependence and in case of no
      timedependence use a polynom with  a_0 = 1 and a_k = 0, k = 1,2,3.
883
    */
884

885 886 887 888
    std::shared_ptr<AbstractTimeDependence> freq_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> ampl_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> phase_atd = nullptr;

889
    dvector_t  unityVec;
890 891 892 893
    unityVec.push_back(1.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
894

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

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

    if (elptr->getPhaseModelName() != "") {
910 911
        phase_atd = AbstractTimeDependence::getTimeDependence(elptr->getPhaseModelName());
        *gmsg << "* Variable phase RF Model name " << elptr->getPhaseModelName() << endl;
912
    }
913 914
    else
        phase_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
915

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

snuverink_j's avatar
snuverink_j committed
919
    double BcParameter[8] = {};
920

921 922
    BcParameter[0] = 0.001 * rmin;
    BcParameter[1] = 0.001 * rmax;
923
    BcParameter[2] = 0.001 * pdis;
gsell's avatar
gsell committed