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>
frey_m's avatar
frey_m committed
25
#include <numeric>
26
#include <cmath>
27 28

#include "AbstractObjects/Element.h"
29
#include "AbstractObjects/OpalData.h"
gsell's avatar
gsell committed
30

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

66 67
#include "Algorithms/Ctunes.h"
#include "Algorithms/PolynomialTimeDependence.h"
gsell's avatar
gsell committed
68 69

#include "Beamlines/Beamline.h"
70
#include "Beamlines/FlaggedBeamline.h"
gsell's avatar
gsell committed
71

72
#include "Elements/OpalBeamline.h"
gsell's avatar
gsell committed
73 74 75 76

#include "Physics/Physics.h"

#include "Utilities/OpalException.h"
77
#include "Utilities/Options.h"
gsell's avatar
gsell committed
78

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

82
#include "Structure/BoundaryGeometry.h"
83
#include "Structure/DataSink.h"
Pedro Calvo Portela's avatar
Pedro Calvo Portela committed
84
#include "Structure/LossDataSink.h"
frey_m's avatar
frey_m committed
85

86

87
constexpr double c_mmtns = Physics::c * 1.0e-6; // m/s --> mm/ns
88

89 90 91
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
92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107

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
108
                                                   PartBunchBase<double, 3> *bunch,
gsell's avatar
gsell committed
109 110 111
                                                   DataSink &ds,
                                                   const PartData &reference,
                                                   bool revBeam, bool revTrack,
frey_m's avatar
frey_m committed
112
                                                   int maxSTEPS, int timeIntegrator,
113 114 115 116 117
                                                   const int& numBunch,
                                                   const double& mbEta,
                                                   const double& mbPara,
                                                   const std::string& mbMode,
                                                   const std::string& mbBinning)
frey_m's avatar
frey_m committed
118 119 120 121 122 123 124 125 126 127
    : Tracker(beamline, bunch, reference, revBeam, revTrack)
    , bgf_m(nullptr)
    , maxSteps_m(maxSTEPS)
    , lastDumpedStep_m(0)
    , myNode_m(Ippl::myNode())
    , initialLocalNum_m(bunch->getLocalNum())
    , initialTotalNum_m(bunch->getTotalNum())
    , opalRing_m(nullptr)
    , itsStepper_mp(nullptr)
{
gsell's avatar
gsell committed
128 129 130
    itsBeamline = dynamic_cast<Beamline *>(beamline.clone());
    itsDataSink = &ds;

131 132 133
    if ( numBunch > 1 ) {
        mbHandler_m = std::unique_ptr<MultiBunchHandler>(
            new MultiBunchHandler(bunch, numBunch, mbEta,
frey_m's avatar
frey_m committed
134
                                  mbPara, mbMode, mbBinning)
135 136 137
        );
    }

gsell's avatar
gsell committed
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
    PluginElemTimer_m  = IpplTimings::getTimer("PluginElements");
    DelParticleTimer_m = IpplTimings::getTimer("DeleteParticles");
144

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

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

/**
 * Destructor ParallelCyclotronTracker
 *
 */
ParallelCyclotronTracker::~ParallelCyclotronTracker() {
Pedro Calvo Portela's avatar
Pedro Calvo Portela committed
170 171
    if(bgf_m)
        lossDs_m->save();
172 173
    for(Component* component : myElements) {
        delete(component);
gsell's avatar
gsell committed
174
    }
175 176
    for(auto fd : FieldDimensions) {
        delete(fd);
gsell's avatar
gsell committed
177 178
    }
    delete itsBeamline;
179
    // delete opalRing_m;
gsell's avatar
gsell committed
180 181
}

frey_m's avatar
frey_m committed
182

183
void ParallelCyclotronTracker::bgf_main_collision_test() {
184
    if(!bgf_m) return;
185

186
    Inform msg("bgf_main_collision_test ");
187

188
    /**
189
     *Here we check if a particle is outside the domain, flag it for deletion
190
     */
191

192
    Vector_t intecoords = 0.0;
193

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

197
    int triId = 0;
frey_m's avatar
frey_m committed
198 199 200
    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);
201
        if(res >= 0) {
ext-calvo_p's avatar
ext-calvo_p committed
202
            lossDs_m->addParticle(itsBunch_m->R[i], itsBunch_m->P[i],
Pedro Calvo Portela's avatar
Pedro Calvo Portela committed
203 204
                                  itsBunch_m->ID[i], itsBunch_m->getT()*1e9,
                                  turnnumber_m, itsBunch_m->bunchNum[i]);
frey_m's avatar
frey_m committed
205
            itsBunch_m->Bin[i] = -1;
206
            Inform gmsgALL("OPAL", INFORM_ALL_NODES);
Pedro Calvo Portela's avatar
Pedro Calvo Portela committed
207
            gmsgALL << level4 << "* Particle " << itsBunch_m->ID[i] << " lost on boundary geometry" << endl;
208
        }
209
    }
210 211
}

212
// only used for dumping into stat file
frey_m's avatar
frey_m committed
213 214
void ParallelCyclotronTracker::dumpAngle(const double& theta,
                                         double& prevAzimuth,
frey_m's avatar
frey_m committed
215 216
                                         double& azimuth,
                                         const short& bunchNr)
frey_m's avatar
frey_m committed
217 218
{
    if ( prevAzimuth < 0.0 ) { // only at first occurrence
frey_m's avatar
frey_m committed
219 220
        double plus = 0.0;
        if ( OpalData::getInstance()->inRestartRun() ) {
frey_m's avatar
frey_m committed
221
            plus = 360.0 * (turnnumber_m - bunchNr);
frey_m's avatar
frey_m committed
222
        }
frey_m's avatar
frey_m committed
223
        azimuth = theta + plus;
224
    } else {
frey_m's avatar
frey_m committed
225
        double dtheta = theta - prevAzimuth;
226 227 228 229 230 231
        if ( dtheta < 0 ) {
            dtheta += 360.0;
        }
        if ( dtheta > 180 ) { // rotating clockwise, reduce angle
            dtheta -= 360;
        }
frey_m's avatar
frey_m committed
232
        azimuth += dtheta;
233
    }
frey_m's avatar
frey_m committed
234
    prevAzimuth = theta;
235 236 237
}


frey_m's avatar
frey_m committed
238 239 240 241 242 243
double ParallelCyclotronTracker::computeRadius(const Vector_t& meanR) const {
    // New OPAL 2.0: m --> mm
    return 1000.0 * std::sqrt(meanR(0) * meanR(0) + meanR(1) * meanR(1));
}


frey_m's avatar
frey_m committed
244 245
void ParallelCyclotronTracker::computePathLengthUpdate(std::vector<double>& dl,
                                                       const double& dt)
246
{
frey_m's avatar
frey_m committed
247 248
    // the last element in dotP is the dot-product over all particles
    std::vector<double> dotP(dl.size());
249
    if ( Options::psDumpFrame == Options::BUNCH_MEAN || isMultiBunch()) {
250

frey_m's avatar
frey_m committed
251 252
        for(unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) {
            dotP[itsBunch_m->bunchNum[i]] += dot(itsBunch_m->P[i], itsBunch_m->P[i]);
253 254
        }

frey_m's avatar
frey_m committed
255
        allreduce(dotP.data(), dotP.size(), std::plus<double>());
256

frey_m's avatar
frey_m committed
257 258 259
        // dot-product over all particles
        double sum = std::accumulate(dotP.begin(), dotP.end() - 1, 0);
        dotP.back() = sum / double(itsBunch_m->getTotalNum());
frey_m's avatar
frey_m committed
260

frey_m's avatar
frey_m committed
261 262 263 264
        // bunch specific --> multi-bunches only
        for (short b = 0; b < (short)dotP.size() - 1; ++b) {
            dotP[b] /= double(itsBunch_m->getTotalNumPerBunch(b));
        }
265 266

    } else if ( itsBunch_m->getLocalNum() == 0 ) {
267
        // here we are in Options::GLOBAL mode
frey_m's avatar
frey_m committed
268
        dotP[0] = 0.0;
269
    } else {
270
        // here we are in Options::GLOBAL mode
frey_m's avatar
frey_m committed
271
        dotP[0] = dot(itsBunch_m->P[0], itsBunch_m->P[0]);
272 273
    }

frey_m's avatar
frey_m committed
274 275 276 277 278
    for (size_t i = 0; i < dotP.size(); ++i) {
        double const gamma = std::sqrt(1.0 + dotP[i]);
        double const beta  = std::sqrt(dotP[i]) / gamma;
        dl[i] = c_mmtns * dt * 1.0e-3 * beta; // unit: m
    }
279 280 281
}


282

gsell's avatar
gsell committed
283 284 285 286 287
/**
 *
 *
 * @param fn Base file name
 */
288
void ParallelCyclotronTracker::openFiles(size_t numFiles, std::string SfileName) {
gsell's avatar
gsell committed
289

290
    for (unsigned int i=0; i<numFiles; i++) {
291 292 293 294 295 296 297 298
        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
299

300 301 302 303 304 305
        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;
306
    }
gsell's avatar
gsell committed
307 308 309 310 311 312 313 314
}

/**
 * Close all files related to
 * special output in the Cyclotron
 * mode.
 */
void ParallelCyclotronTracker::closeFiles() {
315
    for (auto & file : outfTheta_m) {
316
        file->close();
317
    }
gsell's avatar
gsell committed
318 319
}

320
/**
321 322 323
 *
 * @param ring
 */
kraus's avatar
kraus committed
324
void ParallelCyclotronTracker::visitRing(const Ring &ring) {
Daniel Winklehner's avatar
Daniel Winklehner committed
325

326
    *gmsg << "* ----------------------------- Adding Ring ------------------------------ *" << endl;
Daniel Winklehner's avatar
Daniel Winklehner committed
327

snuverink_j's avatar
snuverink_j committed
328
    delete opalRing_m;
Daniel Winklehner's avatar
Daniel Winklehner committed
329

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

332
    myElements.push_back(opalRing_m);
Daniel Winklehner's avatar
Daniel Winklehner committed
333

frey_m's avatar
frey_m committed
334
    opalRing_m->initialise(itsBunch_m);
335 336 337 338

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

340
    if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
kraus's avatar
kraus committed
341
        throw OpalException("Error in ParallelCyclotronTracker::visitRing",
342 343
                            "PHIINIT is out of [-180, 180)!");
    }
Daniel Winklehner's avatar
Daniel Winklehner committed
344 345

    referenceZ = 0.0;
346
    referencePz = 0.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
347 348

    referencePtot = itsReference.getGamma() * itsReference.getBeta();
349
    referencePt = std::sqrt(referencePtot * referencePtot - referencePr * referencePr);
Daniel Winklehner's avatar
Daniel Winklehner committed
350

351 352
    if(referencePtot < 0.0)
        referencePt *= -1.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
353

354 355
    sinRefTheta_m = std::sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = std::cos(referenceTheta * Physics::deg2rad);
gsell's avatar
gsell committed
356

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

359
    buildupFieldList(BcParameter, ElementBase::RING, opalRing_m);
360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375

    // 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
376 377 378 379 380 381 382 383

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

384
    *gmsg << "* -------------------------- Adding Cyclotron ---------------------------- *" << endl;
gsell's avatar
gsell committed
385

386 387
    cycl_m = dynamic_cast<Cyclotron *>(cycl.clone());
    myElements.push_back(cycl_m);
388

389
    // Is this a Spiral Inflector Simulation? If yes, we'll give the user some
390
    // useful information
391
    spiral_flag = cycl_m->getSpiralFlag();
392 393 394 395

    if(spiral_flag) {

        *gmsg << endl << "* This is a Spiral Inflector Simulation! This means the following:" << endl;
396
        *gmsg         << "* 1.) It is up to the user to provide appropriate geometry, electric and magnetic fields!" << endl;
397 398
        *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
399
        *gmsg         << "* 2.) For high currents it is strongly recommended to use the SAAMG fieldsolver," << endl;
400 401
        *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;
402 403
        *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;
404 405 406
        *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;
frey_m's avatar
frey_m committed
407 408 409
        if (isMultiBunch()) {
            mbHandler_m = nullptr;
        }
410 411
    }

412
    // Fresh run (no restart):
413
    if(!OpalData::getInstance()->inRestartRun()) {
414

415
        // Get reference values from cyclotron element
416
        // For now, these are still stored in mm. should be the only ones. -DW
417 418 419 420 421
        referenceR     = cycl_m->getRinit();
        referenceTheta = cycl_m->getPHIinit();
        referenceZ     = cycl_m->getZinit();
        referencePr    = cycl_m->getPRinit();
        referencePz    = cycl_m->getPZinit();
422

423
        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
424 425
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                "PHIINIT is out of [-180, 180)!");
426 427 428 429 430
        }

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

        // Calculate reference azimuthal (tangential) momentum from total-, z- and radial momentum:
431
        float insqrt = referencePtot * referencePtot - \
432
            referencePr * referencePr - referencePz * referencePz;
433 434 435 436 437

        if(insqrt < 0) {

            if(insqrt > -1.0e-10) {

438
                referencePt = 0.0;
439 440 441

            } else {

442
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
443
                                    "Pt imaginary!");
444 445 446 447
            }

        } else {

448
            referencePt = std::sqrt(insqrt);
449 450
        }

451
        if(referencePtot < 0.0)
Daniel Winklehner's avatar
Daniel Winklehner committed
452
            referencePt *= -1.0;
453 454
        // End calculate referencePt

455
        // Restart a run:
456 457
    } else {

458
        // If the user wants to save the restarted run in local frame,
459
        // make sure the previous h5 file was local too
460
      if (Options::psDumpFrame != Options::GLOBAL) {
winklehner_d's avatar
-DW  
winklehner_d committed
461
        if (!previousH5Local) {
462 463
                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
464
        }
465 466
            // 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
467 468
      } else {
        if (previousH5Local) {
469 470
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                    "You are trying a global restart from a local h5 file!");
471
            }
472
        }
473

474
        // Adjust some of the reference variables from the h5 file
475 476
        referencePhi *= Physics::deg2rad;
        referencePsi *= Physics::deg2rad;
477
        referencePtot = bega;
478
        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
479
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
480 481
                                "PHIINIT is out of [-180, 180)!");
        }
482 483
    }

484 485
    sinRefTheta_m = std::sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = std::cos(referenceTheta * Physics::deg2rad);
486

487
    *gmsg << endl;
adelmann's avatar
adelmann committed
488
    *gmsg << "* Bunch global starting position:" << endl;
489 490
    *gmsg << "* RINIT = " << referenceR  << " [mm]" << endl;
    *gmsg << "* PHIINIT = " << referenceTheta << " [deg]" << endl;
491
    *gmsg << "* ZINIT = " << referenceZ << " [mm]" << endl;
492
    *gmsg << endl;
adelmann's avatar
adelmann committed
493
    *gmsg << "* Bunch global starting momenta:" << endl;
gsell's avatar
gsell committed
494 495
    *gmsg << "* Initial gamma = " << itsReference.getGamma() << endl;
    *gmsg << "* Initial beta = " << itsReference.getBeta() << endl;
496
    *gmsg << "* Reference total momentum (beta * gamma) = " << referencePtot * 1000.0 << " [MCU]" << endl;
497 498 499
    *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;
500
    *gmsg << endl;
adelmann's avatar
adelmann committed
501

502
    double sym = cycl_m->getSymmetry();
snuverink_j's avatar
snuverink_j committed
503
    *gmsg << "* " << sym << "-fold field symmetry " << endl;
gsell's avatar
gsell committed
504

505
    // ckr: this just returned the default value as defined in Component.h
506
    // double rff = cycl_m->getRfFrequ(0);
507
    // *gmsg << "* Rf frequency= " << rff << " [MHz]" << endl;
gsell's avatar
gsell committed
508

509
    std::string fmfn = cycl_m->getFieldMapFN();
510
    *gmsg << "* Field map file name = " << fmfn << " " << endl;
gsell's avatar
gsell committed
511

512
    std::string type = cycl_m->getCyclotronType();
513
    *gmsg << "* Type of cyclotron = " << type << " " << endl;
514

515 516
    double rmin = cycl_m->getMinR();
    double rmax = cycl_m->getMaxR();
adelmann's avatar
adelmann committed
517
    *gmsg << "* Radial aperture = " << rmin << " ... " << rmax<<" [m] "<< endl;
518

519 520
    double zmin = cycl_m->getMinZ();
    double zmax = cycl_m->getMaxZ();
adelmann's avatar
adelmann committed
521
    *gmsg << "* Vertical aperture = " << zmin << " ... " << zmax<<" [m]"<< endl;
gsell's avatar
gsell committed
522

523 524
    double h = cycl_m->getCyclHarm();
    *gmsg << "* Number of trimcoils = " << cycl_m->getNumberOfTrimcoils() << endl;
525
    *gmsg << "* Harmonic number h = " << h << " " << endl;
526

527 528 529 530 531 532 533 534 535
    if (type == std::string("BANDRF")) {
        double escale = cycl_m->getEScale(0);
        *gmsg << "* RF field scale factor = " << escale << endl;
        double rfphi= cycl_m->getRfPhi(0);
        *gmsg << "* RF inital phase = " << rfphi * Physics::rad2deg << " [deg]" << endl;
        bool superpose = cycl_m->getSuperpose(0);
        *gmsg << std::boolalpha << "* Superpose electric field maps -> " << superpose << endl;
    }

gsell's avatar
gsell committed
536 537 538 539 540 541 542
    /**
     * 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
543
     * fieldflag = 5, readin FFA format file for MSU/FNAL FFA
gsell's avatar
gsell committed
544
     * fieldflag = 6, readin both median plane B field map and 3D E field map of RF cavity for compact cyclotron
545
     * fieldflag = 7, read in fields for Daniel's synchrocyclotron simulations
gsell's avatar
gsell committed
546
     */
547
    int fieldflag = cycl_m->getFieldFlag(type);
gsell's avatar
gsell committed
548

549
    // Read in cyclotron field maps (midplane + 3D fields if desired).
550
    cycl_m->initialise(itsBunch_m, fieldflag, cycl_m->getBScale());
gsell's avatar
gsell committed
551

snuverink_j's avatar
snuverink_j committed
552
    double BcParameter[8] = {};
553

554 555
    BcParameter[0] = 0.001 * cycl_m->getRmin();
    BcParameter[1] = 0.001 * cycl_m->getRmax();
gsell's avatar
gsell committed
556

557
    // Store inner radius and outer radius of cyclotron field map in the list
558
    buildupFieldList(BcParameter, ElementBase::CYCLOTRON, cycl_m);
gsell's avatar
gsell committed
559 560 561 562 563 564 565 566 567 568
}

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

569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601

void ParallelCyclotronTracker::visitBeamStripping(const BeamStripping &bstp) {
    *gmsg << "* ------------------------------ Beam Stripping ------------------------------" << endl;

    BeamStripping* elptr = dynamic_cast<BeamStripping *>(bstp.clone());
    myElements.push_back(elptr);

    double BcParameter[8] = {};

    if(elptr->getPressureMapFN() == "") {
        double pressure = elptr->getPressure();
        *gmsg << "* Pressure     = " << pressure << " [mbar]" << endl;
        BcParameter[0] = pressure;
    }

    double temperature = elptr->getTemperature();
    *gmsg << "* Temperature  = " << temperature << " [K]" << endl;

    std::string gas = elptr->getResidualGas();
    *gmsg << "* Residual gas = " << gas << endl;

    bool stop = elptr->getStop();
    *gmsg << std::boolalpha << "* Particles stripped will be deleted after interaction -> " << stop << endl;

    elptr->initialise(itsBunch_m, elptr->getPScale());
    
    BcParameter[1] = temperature;
    BcParameter[2] = stop;

    buildupFieldList(BcParameter, ElementBase::BEAMSTRIPPING, elptr);

}

gsell's avatar
gsell committed
602 603 604 605 606
/**
 *
 *
 * @param coll
 */
607
void ParallelCyclotronTracker::visitCCollimator(const CCollimator &coll) {
gsell's avatar
gsell committed
608

609
    *gmsg << "* ------------------------------ Collimator ------------------------------" << endl;
gsell's avatar
gsell committed
610

611
    CCollimator* elptr = dynamic_cast<CCollimator *>(coll.clone());
612
    myElements.push_back(elptr);
gsell's avatar
gsell committed
613

614
    double xstart = elptr->getXStart();
615
    *gmsg << "* Xstart  = " << xstart << " [mm]" << endl;
gsell's avatar
gsell committed
616

617
    double xend = elptr->getXEnd();
618
    *gmsg << "* Xend    = " << xend << " [mm]" << endl;
gsell's avatar
gsell committed
619

620
    double ystart = elptr->getYStart();
621
    *gmsg << "* Ystart  = " << ystart << " [mm]" << endl;
gsell's avatar
gsell committed
622

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

626
    double zstart = elptr->getZStart();
627
    *gmsg << "* Zstart  = " << zstart << " [mm]" << endl;
628 629

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

632
    double width = elptr->getWidth();
633
    *gmsg << "* Width   = " << width << " [mm]" << endl;
gsell's avatar
gsell committed
634

frey_m's avatar
frey_m committed
635
    elptr->initialise(itsBunch_m);
gsell's avatar
gsell committed
636

snuverink_j's avatar
snuverink_j committed
637
    double BcParameter[8] = {};
638

639 640 641 642 643 644
    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 ;

645
    buildupFieldList(BcParameter, ElementBase::CCOLLIMATOR, elptr);
gsell's avatar
gsell committed
646 647 648 649 650 651 652 653 654 655 656 657
}

/**
 *
 *
 * @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
658 659 660 661 662 663 664 665 666 667 668 669
/**
 *
 *
 * @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
670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689
/**
 *
 *
 * @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()));
}

690 691
/**
 *
692 693
 *
 *  @param
694 695 696 697 698
 */
void ParallelCyclotronTracker::visitFlexibleCollimator(const FlexibleCollimator &) {

}

gsell's avatar
gsell committed
699 700 701 702 703 704 705 706 707 708
/**
 *
 *
 * @param lamb
 */
void ParallelCyclotronTracker::visitLambertson(const Lambertson &lamb) {
    *gmsg << "In Lambertson; L= " << lamb.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Lambertson *>(lamb.clone()));
}

709 710 711
void ParallelCyclotronTracker::visitOffset(const Offset & off) {
    if (opalRing_m == NULL)
        throw OpalException(
712 713
                            "ParallelCylcotronTracker::visitOffset",
                            "Attempt to place an offset when Ring not defined");
714 715
    Offset* offNonConst = const_cast<Offset*>(&off);
    offNonConst->updateGeometry(opalRing_m->getNextPosition(),
716
                                opalRing_m->getNextNormal());
717 718 719
    opalRing_m->appendElement(off);
}

gsell's avatar
gsell committed
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
/**
 *
 *
 * @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
752 753 754 755 756 757 758
/**
 *
 *
 * @param multT
 */
void ParallelCyclotronTracker::visitMultipoleT(const MultipoleT &multT) {
    *gmsg << "Adding MultipoleT" << endl;
759
    if (opalRing_m != NULL) {
ext-rogers_c's avatar
ext-rogers_c committed
760
        opalRing_m->appendElement(multT);
761
    } else {
ext-rogers_c's avatar
ext-rogers_c committed
762 763
        throw OpalException("ParallelCyclotronTracker::visitMultipoleT",
                            "Need to define a RINGDEFINITION to use MultipoleT element");
764
    }
ext-rogers_c's avatar
ext-rogers_c committed
765 766 767
    myElements.push_back(dynamic_cast<MultipoleT *>(multT.clone()));
}

768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815
/**
 *
 *
 * @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
816 817 818 819 820 821
/**
 *
 *
 * @param prob
 */
void ParallelCyclotronTracker::visitProbe(const Probe &prob) {
822
    *gmsg << "* ------------------------------  Probe ------------------------------" << endl;
823 824
    Probe *elptr = dynamic_cast<Probe *>(prob.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
825

826 827
    *gmsg << "* Name    = " << elptr->getName() << endl;

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

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

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

837
    double yend = elptr->getYEnd();
838
    *gmsg << "* YEnd    = " << yend << " [mm]" << endl;
gsell's avatar
gsell committed
839 840

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

snuverink_j's avatar
snuverink_j committed
843
    double BcParameter[8] = {};
844

845 846 847 848
    BcParameter[0] = 0.001 * xstart ;
    BcParameter[1] = 0.001 * xend;
    BcParameter[2] = 0.001 * ystart ;
    BcParameter[3] = 0.001 * yend;
849
    BcParameter[4] = 0.001 ; // width
gsell's avatar
gsell committed
850 851

    // store probe parameters in the list
852
    buildupFieldList(BcParameter, ElementBase::PROBE, elptr);
gsell's avatar
gsell committed
853 854 855 856 857 858 859 860 861 862 863 864
}

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

865 866 867 868 869 870
void ParallelCyclotronTracker::visitSBend3D(const SBend3D &bend) {
    *gmsg << "Adding SBend3D" << endl;
    if (opalRing_m != NULL)
        opalRing_m->appendElement(bend);
    else
        throw OpalException("ParallelCyclotronTracker::visitSBend3D",
871
                            "Need to define a RINGDEFINITION to use SBend3D element");
872 873
}

ext-rogers_c's avatar
ext-rogers_c committed
874 875
void ParallelCyclotronTracker::visitScalingFFAMagnet(const ScalingFFAMagnet &bend) {
    *gmsg << "Adding ScalingFFAMagnet" << endl;
ext-rogers_c's avatar
ext-rogers_c committed
876
    if (opalRing_m != NULL) {
877
        opalRing_m->appendElement(bend);
ext-rogers_c's avatar
ext-rogers_c committed
878
    } else {
ext-rogers_c's avatar
ext-rogers_c committed
879 880
        throw OpalException("ParallelCyclotronTracker::visitScalingFFAMagnet",
                            "Need to define a RINGDEFINITION to use ScalingFFAMagnet element");
ext-rogers_c's avatar
ext-rogers_c committed
881
    }
882 883
}

884 885 886 887 888 889
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",
890
                            "Need to define a RINGDEFINITION to use VariableRFCavity element");
891 892 893 894 895 896 897 898 899 900
}

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

903 904 905 906 907 908 909 910 911
void ParallelCyclotronTracker::visitVerticalFFAMagnet(