ParallelCyclotronTracker.cpp 126 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 123
    : Tracker(beamline, bunch, reference, revBeam, revTrack)
    , itsMBDump_m(new MultiBunchDump())
    , bgf_m(nullptr)
    , maxSteps_m(maxSTEPS)
124
    , numBunch_m(numBunch)
frey_m's avatar
frey_m committed
125 126 127 128 129 130 131 132 133
    , 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
134 135
    itsBeamline = dynamic_cast<Beamline *>(beamline.clone());
    itsDataSink = &ds;
136
    multiBunchMode_m = MB_MODE::NONE;
gsell's avatar
gsell committed
137 138 139 140 141

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

143 144 145 146 147 148 149 150 151
    // 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;
152

153 154 155 156 157 158 159 160
    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
161 162 163 164 165 166 167
}

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

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

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

frey_m's avatar
frey_m committed
202 203 204 205 206 207 208 209 210 211 212 213
    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.");
    }
}

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

222
    Inform msg("bgf_main_collision_test ");
223

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

228
    Vector_t intecoords = 0.0;
229

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

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


gsell's avatar
gsell committed
244 245 246 247 248
/**
 *
 *
 * @param fn Base file name
 */
249
void ParallelCyclotronTracker::openFiles(size_t numFiles, std::string SfileName) {
gsell's avatar
gsell committed
250

251
    for (unsigned int i=0; i<numFiles; i++) {
252 253 254 255 256 257 258 259
        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
260

261 262 263 264 265 266
        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;
267
    }
gsell's avatar
gsell committed
268 269 270 271 272 273 274 275
}

/**
 * Close all files related to
 * special output in the Cyclotron
 * mode.
 */
void ParallelCyclotronTracker::closeFiles() {
276
    for (auto & file : outfTheta_m) {
277
        file->close();
278
    }
gsell's avatar
gsell committed
279 280
}

281
/**
282 283 284
 *
 * @param ring
 */
kraus's avatar
kraus committed
285
void ParallelCyclotronTracker::visitRing(const Ring &ring) {
Daniel Winklehner's avatar
Daniel Winklehner committed
286

287
    *gmsg << "* ----------------------------- Adding Ring ------------------------------ *" << endl;
Daniel Winklehner's avatar
Daniel Winklehner committed
288

snuverink_j's avatar
snuverink_j committed
289
    delete opalRing_m;
Daniel Winklehner's avatar
Daniel Winklehner committed
290

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

293
    myElements.push_back(opalRing_m);
Daniel Winklehner's avatar
Daniel Winklehner committed
294

frey_m's avatar
frey_m committed
295
    opalRing_m->initialise(itsBunch_m);
296 297 298 299

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

301
    if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
kraus's avatar
kraus committed
302
        throw OpalException("Error in ParallelCyclotronTracker::visitRing",
303 304
                            "PHIINIT is out of [-180, 180)!");
    }
Daniel Winklehner's avatar
Daniel Winklehner committed
305 306

    referenceZ = 0.0;
307
    referencePz = 0.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
308 309 310 311

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

312 313
    if(referencePtot < 0.0)
        referencePt *= -1.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
314

315 316
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
gsell's avatar
gsell committed
317

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

320
    buildupFieldList(BcParameter, ElementBase::RING, opalRing_m);
321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336

    // 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
337 338 339 340 341 342 343 344

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

345
    *gmsg << "* -------------------------- Adding Cyclotron ---------------------------- *" << endl;
gsell's avatar
gsell committed
346

347 348
    Cyclotron *elptr = dynamic_cast<Cyclotron *>(cycl.clone());
    myElements.push_back(elptr);
349

350
    // Is this a Spiral Inflector Simulation? If yes, we'll give the user some
351 352 353 354 355 356
    // useful information
    spiral_flag = elptr->getSpiralFlag();

    if(spiral_flag) {

        *gmsg << endl << "* This is a Spiral Inflector Simulation! This means the following:" << endl;
357
        *gmsg         << "* 1.) It is up to the user to provide appropriate geometry, electric and magnetic fields!" << endl;
358 359
        *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
360
        *gmsg         << "* 2.) For high currents it is strongly recommended to use the SAAMG fieldsolver," << endl;
361 362
        *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;
363 364
        *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;
365 366 367 368 369 370 371
        *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;

    }

372
    // Fresh run (no restart):
373
    if(!OpalData::getInstance()->inRestartRun()) {
374

375
        // Get reference values from cyclotron element
376
        // For now, these are still stored in mm. should be the only ones. -DW
377 378 379 380
        referenceR     = elptr->getRinit();
        referenceTheta = elptr->getPHIinit();
        referenceZ     = elptr->getZinit();
        referencePr    = elptr->getPRinit();
381 382
        referencePz    = elptr->getPZinit();

383
        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
384 385
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                "PHIINIT is out of [-180, 180)!");
386 387 388 389 390
        }

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

        // Calculate reference azimuthal (tangential) momentum from total-, z- and radial momentum:
391
        float insqrt = referencePtot * referencePtot - \
392
            referencePr * referencePr - referencePz * referencePz;
393 394 395 396 397

        if(insqrt < 0) {

            if(insqrt > -1.0e-10) {

398
                referencePt = 0.0;
399 400 401

            } else {

402
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
403
                                    "Pt imaginary!");
404 405 406 407 408 409 410
            }

        } else {

            referencePt = sqrt(insqrt);
        }

411
        if(referencePtot < 0.0)
Daniel Winklehner's avatar
Daniel Winklehner committed
412
            referencePt *= -1.0;
413 414
        // End calculate referencePt

415
        // Restart a run:
416 417
    } else {

418
        // If the user wants to save the restarted run in local frame,
419
        // make sure the previous h5 file was local too
420
      if (Options::psDumpFrame != Options::GLOBAL) {
winklehner_d's avatar
-DW  
winklehner_d committed
421
        if (!previousH5Local) {
422 423
                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
424
        }
425 426
            // 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
427 428
      } else {
        if (previousH5Local) {
429 430
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                    "You are trying a global restart from a local h5 file!");
431
            }
432
        }
433

434
        // Adjust some of the reference variables from the h5 file
435 436
        referencePhi *= Physics::deg2rad;
        referencePsi *= Physics::deg2rad;
437
        referencePtot = bega;
438
        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
439
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
440 441
                                "PHIINIT is out of [-180, 180)!");
        }
442 443
    }

444 445
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
446

447
    *gmsg << endl;
adelmann's avatar
adelmann committed
448
    *gmsg << "* Bunch global starting position:" << endl;
449 450
    *gmsg << "* RINIT = " << referenceR  << " [mm]" << endl;
    *gmsg << "* PHIINIT = " << referenceTheta << " [deg]" << endl;
451
    *gmsg << "* ZINIT = " << referenceZ << " [mm]" << endl;
452
    *gmsg << endl;
adelmann's avatar
adelmann committed
453
    *gmsg << "* Bunch global starting momenta:" << endl;
gsell's avatar
gsell committed
454 455
    *gmsg << "* Initial gamma = " << itsReference.getGamma() << endl;
    *gmsg << "* Initial beta = " << itsReference.getBeta() << endl;
456
    *gmsg << "* Reference total momentum (beta * gamma) = " << referencePtot * 1000.0 << " [MCU]" << endl;
457 458 459
    *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;
460
    *gmsg << endl;
adelmann's avatar
adelmann committed
461

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

465 466 467
    // 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
468

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

472
    std::string type = elptr->getCyclotronType();
473
    *gmsg << "* Type of cyclotron = " << type << " " << endl;
474

475 476
    double rmin = elptr->getMinR();
    double rmax = elptr->getMaxR();
adelmann's avatar
adelmann committed
477
    *gmsg << "* Radial aperture = " << rmin << " ... " << rmax<<" [m] "<< endl;
478 479 480

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

    double h = elptr->getCyclHarm();
484
    *gmsg << "* Number of trimcoils = " << elptr->getNumberOfTrimcoils() << endl;
485
    *gmsg << "* Harmonic number h = " << h << " " << endl;
486

gsell's avatar
gsell committed
487 488 489 490 491 492 493
    /**
     * 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
494
     * fieldflag = 5, readin FFA format file for MSU/FNAL FFA
gsell's avatar
gsell committed
495
     * fieldflag = 6, readin both median plane B field map and 3D E field map of RF cavity for compact cyclotron
496
     * fieldflag = 7, read in fields for Daniel's synchrocyclotron simulations
gsell's avatar
gsell committed
497
     */
498
    int fieldflag = elptr->getFieldFlag(type);
gsell's avatar
gsell committed
499

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

snuverink_j's avatar
snuverink_j committed
503
    double BcParameter[8] = {};
504

505 506
    BcParameter[0] = 0.001 * elptr->getRmin();
    BcParameter[1] = 0.001 * elptr->getRmax();
gsell's avatar
gsell committed
507

508
    // Store inner radius and outer radius of cyclotron field map in the list
509
    buildupFieldList(BcParameter, ElementBase::CYCLOTRON, elptr);
gsell's avatar
gsell committed
510 511 512 513 514 515 516 517 518 519 520 521 522 523 524
}

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

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

527
    *gmsg << "* --------- Collimator -----------------------------" << endl;
gsell's avatar
gsell committed
528

529
    CCollimator* elptr = dynamic_cast<CCollimator *>(coll.clone());
530
    myElements.push_back(elptr);
gsell's avatar
gsell committed
531

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

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

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

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

544
    double zstart = elptr->getZStart();
adelmann's avatar
adelmann committed
545
    *gmsg << "* Zstart= " << zstart << " [mm]" << endl;
546 547

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

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

frey_m's avatar
frey_m committed
553
    elptr->initialise(itsBunch_m);
gsell's avatar
gsell committed
554

snuverink_j's avatar
snuverink_j committed
555
    double BcParameter[8] = {};
556

557 558 559 560 561 562
    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 ;

563
    buildupFieldList(BcParameter, ElementBase::CCOLLIMATOR, elptr);
gsell's avatar
gsell committed
564 565 566 567 568 569 570 571 572 573 574 575
}

/**
 *
 *
 * @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
576 577 578 579 580 581 582 583 584 585 586 587
/**
 *
 *
 * @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
588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607
/**
 *
 *
 * @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()));
}

608 609
/**
 *
610 611
 *
 *  @param
612 613 614 615 616
 */
void ParallelCyclotronTracker::visitFlexibleCollimator(const FlexibleCollimator &) {

}

gsell's avatar
gsell committed
617 618 619 620 621 622 623 624 625 626
/**
 *
 *
 * @param lamb
 */
void ParallelCyclotronTracker::visitLambertson(const Lambertson &lamb) {
    *gmsg << "In Lambertson; L= " << lamb.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Lambertson *>(lamb.clone()));
}

627 628 629
void ParallelCyclotronTracker::visitOffset(const Offset & off) {
    if (opalRing_m == NULL)
        throw OpalException(
630 631
                            "ParallelCylcotronTracker::visitOffset",
                            "Attempt to place an offset when Ring not defined");
632 633
    Offset* offNonConst = const_cast<Offset*>(&off);
    offNonConst->updateGeometry(opalRing_m->getNextPosition(),
634
                                opalRing_m->getNextNormal());
635 636 637
    opalRing_m->appendElement(off);
}

gsell's avatar
gsell committed
638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669
/**
 *
 *
 * @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
670 671 672 673 674 675 676
/**
 *
 *
 * @param multT
 */
void ParallelCyclotronTracker::visitMultipoleT(const MultipoleT &multT) {
    *gmsg << "Adding MultipoleT" << endl;
677
    if (opalRing_m != NULL) {
ext-rogers_c's avatar
ext-rogers_c committed
678
        opalRing_m->appendElement(multT);
679
    } else {
ext-rogers_c's avatar
ext-rogers_c committed
680 681
        throw OpalException("ParallelCyclotronTracker::visitMultipoleT",
                            "Need to define a RINGDEFINITION to use MultipoleT element");
682
    }
ext-rogers_c's avatar
ext-rogers_c committed
683 684 685
    myElements.push_back(dynamic_cast<MultipoleT *>(multT.clone()));
}

686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733
/**
 *
 *
 * @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
734 735 736 737 738 739
/**
 *
 *
 * @param prob
 */
void ParallelCyclotronTracker::visitProbe(const Probe &prob) {
740
    *gmsg << "* -----------  Probe -------------------------------" << endl;
741 742
    Probe *elptr = dynamic_cast<Probe *>(prob.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
743

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

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

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

753
    double yend = elptr->getYEnd();
754
    *gmsg << "YEnd= " << yend << " [mm]" << endl;
gsell's avatar
gsell committed
755 756

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

snuverink_j's avatar
snuverink_j committed
759
    double BcParameter[8] = {};
760

761 762 763 764
    BcParameter[0] = 0.001 * xstart ;
    BcParameter[1] = 0.001 * xend;
    BcParameter[2] = 0.001 * ystart ;
    BcParameter[3] = 0.001 * yend;
765
    BcParameter[4] = 0.001 ; // width
gsell's avatar
gsell committed
766 767

    // store probe parameters in the list
768
    buildupFieldList(BcParameter, ElementBase::PROBE, elptr);
gsell's avatar
gsell committed
769 770 771 772 773 774 775 776 777 778 779 780
}

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

781 782 783 784 785 786
void ParallelCyclotronTracker::visitSBend3D(const SBend3D &bend) {
    *gmsg << "Adding SBend3D" << endl;
    if (opalRing_m != NULL)
        opalRing_m->appendElement(bend);
    else
        throw OpalException("ParallelCyclotronTracker::visitSBend3D",
787
                            "Need to define a RINGDEFINITION to use SBend3D element");
788 789
}

ext-rogers_c's avatar
ext-rogers_c committed
790 791
void ParallelCyclotronTracker::visitScalingFFAMagnet(const ScalingFFAMagnet &bend) {
    *gmsg << "Adding ScalingFFAMagnet" << endl;
ext-rogers_c's avatar
ext-rogers_c committed
792
    if (opalRing_m != NULL) {
793
        opalRing_m->appendElement(bend);
ext-rogers_c's avatar
ext-rogers_c committed
794
    } else {
ext-rogers_c's avatar
ext-rogers_c committed
795 796
        throw OpalException("ParallelCyclotronTracker::visitScalingFFAMagnet",
                            "Need to define a RINGDEFINITION to use ScalingFFAMagnet element");
ext-rogers_c's avatar
ext-rogers_c committed
797
    }
798 799
}

800 801 802 803 804 805
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",
806
                            "Need to define a RINGDEFINITION to use VariableRFCavity element");
807 808 809 810 811 812 813 814 815 816
}

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

gsell's avatar
gsell committed
819 820 821 822 823 824 825 826
/**
 *
 *
 * @param as
 */
void ParallelCyclotronTracker::visitRFCavity(const RFCavity &as) {

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

828 829
    RFCavity *elptr = dynamic_cast<RFCavity *>(as.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
830 831 832 833 834 835 836

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

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

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

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

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

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

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

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

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

861 862 863
    /*
      Setup time dependence and in case of no
      timedependence use a polynom with  a_0 = 1 and a_k = 0, k = 1,2,3.
864
    */
865

866 867 868 869
    std::shared_ptr<AbstractTimeDependence> freq_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> ampl_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> phase_atd = nullptr;

870
    dvector_t  unityVec;
871 872 873 874
    unityVec.push_back(1.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
875

876
    if (elptr->getFrequencyModelName() != "") {
877 878
        freq_atd = AbstractTimeDependence::getTimeDependence(elptr->getFrequencyModelName());
        *gmsg << "* Variable frequency RF Model name " << elptr->getFrequencyModelName() << endl;
879
    }
880 881
    else
        freq_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
882 883

    if (elptr->getAmplitudeModelName() != "") {
884 885
        ampl_atd = AbstractTimeDependence::getTimeDependence(elptr->getAmplitudeModelName());
        *gmsg << "* Variable amplitude RF Model name " << elptr->getAmplitudeModelName() << endl;
886
    }
887 888
    else
        ampl_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
889 890

    if (elptr->getPhaseModelName() != "") {
891 892
        phase_atd = AbstractTimeDependence::getTimeDependence(elptr->getPhaseModelName());
        *gmsg << "* Variable phase RF Model name " << elptr->getPhaseModelName() << endl;
893
    }
894 895
    else
        phase_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
896

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

snuverink_j's avatar
snuverink_j committed
900
    double BcParameter[8] = {};
901

902 903
    BcParameter[0] = 0.001 * rmin;
    BcParameter[1] = 0.001 * rmax;
904
    BcParameter[2] = 0.001 * pdis;
gsell's avatar
gsell committed
905 906
    BcParameter[3] = angle;

907
    buildupFieldList(BcParameter, ElementBase::RFCAVITY, elptr);
gsell's avatar
gsell committed
908 909 910 911 912 913 914 915
}

/**
 *
 *
 * @param rfq
 */
void ParallelCyclotronTracker::visitRFQuadrupole(const RFQuadrupole &rfq) {
916
    *gmsg << "In RFQuadrupole; L = " << rfq.getElementLength() << " however the element is missing " << endl;
gsell's avatar
gsell committed
917 918 919 920 921 922 923 924 925
    myElements.push_back(dynamic_cast<RFQuadrupole *>(