ParallelCyclotronTracker.cpp 127 KB
Newer Older
gsell's avatar
gsell committed
1
//
2 3
// Class ParallelCyclotronTracker
//   Tracker for OPAL-Cycl
gsell's avatar
gsell committed
4
//
5 6 7 8 9
// Copyright (c) 2007 - 2014, Jianjun Yang, Andreas Adelmann and Matthias Toggweiler,
//                            Paul Scherrer Institut, Villigen PSI, Switzerland
// Copyright (c) 2014,        Daniel Winklehner, MIT, Cambridge, MA, USA
// Copyright (c) 2012 - 2020, Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved
gsell's avatar
gsell committed
10
//
11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
// Implemented as part of the PhD thesis
// "Beam dynamics in high intensity cyclotrons including neighboring bunch effects"
// and the paper
// "Beam dynamics in high intensity cyclotrons including neighboring bunch effects:
// Model, implementation, and application"
// (https://journals.aps.org/prab/pdf/10.1103/PhysRevSTAB.13.064201)
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
gsell's avatar
gsell committed
27
//
kraus's avatar
kraus committed
28
#include "Algorithms/ParallelCyclotronTracker.h"
29

gsell's avatar
gsell committed
30
#include <fstream>
31 32
#include <iostream>
#include <limits>
gsell's avatar
gsell committed
33
#include <vector>
frey_m's avatar
frey_m committed
34
#include <numeric>
35
#include <cmath>
36 37

#include "AbstractObjects/Element.h"
38
#include "AbstractObjects/OpalData.h"
gsell's avatar
gsell committed
39

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

75 76
#include "Algorithms/Ctunes.h"
#include "Algorithms/PolynomialTimeDependence.h"
gsell's avatar
gsell committed
77 78

#include "Beamlines/Beamline.h"
79
#include "Beamlines/FlaggedBeamline.h"
gsell's avatar
gsell committed
80

81
#include "Elements/OpalBeamline.h"
gsell's avatar
gsell committed
82 83 84 85

#include "Physics/Physics.h"

#include "Utilities/OpalException.h"
86
#include "Utilities/Options.h"
gsell's avatar
gsell committed
87

88
#include "BasicActions/DumpFields.h"
89
#include "BasicActions/DumpEMFields.h"
90

91
#include "Structure/BoundaryGeometry.h"
92
#include "Structure/DataSink.h"
Pedro Calvo Portela's avatar
Pedro Calvo Portela committed
93
#include "Structure/LossDataSink.h"
frey_m's avatar
frey_m committed
94

95

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

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

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
117
                                                   PartBunchBase<double, 3> *bunch,
gsell's avatar
gsell committed
118 119 120
                                                   DataSink &ds,
                                                   const PartData &reference,
                                                   bool revBeam, bool revTrack,
frey_m's avatar
frey_m committed
121
                                                   int maxSTEPS, int timeIntegrator,
122 123 124 125 126
                                                   const int& numBunch,
                                                   const double& mbEta,
                                                   const double& mbPara,
                                                   const std::string& mbMode,
                                                   const std::string& mbBinning)
frey_m's avatar
frey_m committed
127 128 129 130 131 132 133 134 135 136
    : 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
137 138 139
    itsBeamline = dynamic_cast<Beamline *>(beamline.clone());
    itsDataSink = &ds;

140 141 142
    if ( numBunch > 1 ) {
        mbHandler_m = std::unique_ptr<MultiBunchHandler>(
            new MultiBunchHandler(bunch, numBunch, mbEta,
frey_m's avatar
frey_m committed
143
                                  mbPara, mbMode, mbBinning)
144 145 146
        );
    }

gsell's avatar
gsell committed
147 148 149 150
    IntegrationTimer_m = IpplTimings::getTimer("Integration");
    TransformTimer_m   = IpplTimings::getTimer("Frametransform");
    DumpTimer_m        = IpplTimings::getTimer("Dump");
    BinRepartTimer_m   = IpplTimings::getTimer("Binaryrepart");
151 152
    PluginElemTimer_m  = IpplTimings::getTimer("PluginElements");
    DelParticleTimer_m = IpplTimings::getTimer("DeleteParticles");
153

154 155 156 157 158 159 160 161 162
    // 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;
163

164 165 166 167 168 169 170 171
    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
172 173 174 175 176 177 178
}

/**
 * Destructor ParallelCyclotronTracker
 *
 */
ParallelCyclotronTracker::~ParallelCyclotronTracker() {
Pedro Calvo Portela's avatar
Pedro Calvo Portela committed
179 180
    if(bgf_m)
        lossDs_m->save();
181 182
    for(Component* component : myElements) {
        delete(component);
gsell's avatar
gsell committed
183
    }
184 185
    for(auto fd : FieldDimensions) {
        delete(fd);
gsell's avatar
gsell committed
186 187
    }
    delete itsBeamline;
188
    // delete opalRing_m;
gsell's avatar
gsell committed
189 190
}

frey_m's avatar
frey_m committed
191

192
void ParallelCyclotronTracker::bgf_main_collision_test() {
193
    if(!bgf_m) return;
194

195
    Inform msg("bgf_main_collision_test ");
196

197
    /**
198
     *Here we check if a particle is outside the domain, flag it for deletion
199
     */
200

201
    Vector_t intecoords = 0.0;
202

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

206
    int triId = 0;
frey_m's avatar
frey_m committed
207 208 209
    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);
210
        if(res >= 0) {
ext-calvo_p's avatar
ext-calvo_p committed
211
            lossDs_m->addParticle(itsBunch_m->R[i], itsBunch_m->P[i],
Pedro Calvo Portela's avatar
Pedro Calvo Portela committed
212 213
                                  itsBunch_m->ID[i], itsBunch_m->getT()*1e9,
                                  turnnumber_m, itsBunch_m->bunchNum[i]);
frey_m's avatar
frey_m committed
214
            itsBunch_m->Bin[i] = -1;
215
            Inform gmsgALL("OPAL", INFORM_ALL_NODES);
Pedro Calvo Portela's avatar
Pedro Calvo Portela committed
216
            gmsgALL << level4 << "* Particle " << itsBunch_m->ID[i] << " lost on boundary geometry" << endl;
217
        }
218
    }
219 220
}

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


frey_m's avatar
frey_m committed
247 248 249 250 251 252
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
253 254
void ParallelCyclotronTracker::computePathLengthUpdate(std::vector<double>& dl,
                                                       const double& dt)
255
{
frey_m's avatar
frey_m committed
256 257
    // the last element in dotP is the dot-product over all particles
    std::vector<double> dotP(dl.size());
258
    if ( Options::psDumpFrame == Options::BUNCH_MEAN || isMultiBunch()) {
259

frey_m's avatar
frey_m committed
260 261
        for(unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) {
            dotP[itsBunch_m->bunchNum[i]] += dot(itsBunch_m->P[i], itsBunch_m->P[i]);
262 263
        }

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

frey_m's avatar
frey_m committed
266 267 268
        // 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
269

frey_m's avatar
frey_m committed
270 271 272 273
        // bunch specific --> multi-bunches only
        for (short b = 0; b < (short)dotP.size() - 1; ++b) {
            dotP[b] /= double(itsBunch_m->getTotalNumPerBunch(b));
        }
274 275

    } else if ( itsBunch_m->getLocalNum() == 0 ) {
276
        // here we are in Options::GLOBAL mode
frey_m's avatar
frey_m committed
277
        dotP[0] = 0.0;
278
    } else {
279
        // here we are in Options::GLOBAL mode
frey_m's avatar
frey_m committed
280
        dotP[0] = dot(itsBunch_m->P[0], itsBunch_m->P[0]);
281 282
    }

frey_m's avatar
frey_m committed
283 284 285 286 287
    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
    }
288 289 290
}


291

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

299
    for (unsigned int i=0; i<numFiles; i++) {
300 301 302 303 304 305 306 307
        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
308

309 310 311 312 313 314
        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;
315
    }
gsell's avatar
gsell committed
316 317 318 319 320 321 322 323
}

/**
 * Close all files related to
 * special output in the Cyclotron
 * mode.
 */
void ParallelCyclotronTracker::closeFiles() {
324
    for (auto & file : outfTheta_m) {
325
        file->close();
326
    }
gsell's avatar
gsell committed
327 328
}

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

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

snuverink_j's avatar
snuverink_j committed
337
    delete opalRing_m;
Daniel Winklehner's avatar
Daniel Winklehner committed
338

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

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

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

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

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

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

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

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

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

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

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

    // 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
385 386 387 388 389 390 391 392

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

393
    *gmsg << "* -------------------------- Adding Cyclotron ---------------------------- *" << endl;
gsell's avatar
gsell committed
394

395 396
    cycl_m = dynamic_cast<Cyclotron *>(cycl.clone());
    myElements.push_back(cycl_m);
397

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

    if(spiral_flag) {

        *gmsg << endl << "* This is a Spiral Inflector Simulation! This means the following:" << endl;
405
        *gmsg         << "* 1.) It is up to the user to provide appropriate geometry, electric and magnetic fields!" << endl;
406 407
        *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
408
        *gmsg         << "* 2.) For high currents it is strongly recommended to use the SAAMG fieldsolver," << endl;
409 410
        *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;
411 412
        *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;
413 414 415
        *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
416 417 418
        if (isMultiBunch()) {
            mbHandler_m = nullptr;
        }
419 420
    }

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

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

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

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

        // Calculate reference azimuthal (tangential) momentum from total-, z- and radial momentum:
440
        float insqrt = referencePtot * referencePtot - \
441
            referencePr * referencePr - referencePz * referencePz;
442 443 444 445 446

        if(insqrt < 0) {

            if(insqrt > -1.0e-10) {

447
                referencePt = 0.0;
448 449 450

            } else {

451
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
452
                                    "Pt imaginary!");
453 454 455 456
            }

        } else {

457
            referencePt = std::sqrt(insqrt);
458 459
        }

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

464
        // Restart a run:
465 466
    } else {

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

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

493 494
    sinRefTheta_m = std::sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = std::cos(referenceTheta * Physics::deg2rad);
495

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

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

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

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

521
    std::string type = cycl_m->getCyclotronType();
522
    *gmsg << "* Type of cyclotron = " << type << " " << endl;
523

524 525
    double rmin = cycl_m->getMinR();
    double rmax = cycl_m->getMaxR();
adelmann's avatar
adelmann committed
526
    *gmsg << "* Radial aperture = " << rmin << " ... " << rmax<<" [m] "<< endl;
527

528 529
    double zmin = cycl_m->getMinZ();
    double zmax = cycl_m->getMaxZ();
adelmann's avatar
adelmann committed
530
    *gmsg << "* Vertical aperture = " << zmin << " ... " << zmax<<" [m]"<< endl;
gsell's avatar
gsell committed
531

532 533
    double h = cycl_m->getCyclHarm();
    *gmsg << "* Number of trimcoils = " << cycl_m->getNumberOfTrimcoils() << endl;
534
    *gmsg << "* Harmonic number h = " << h << " " << endl;
535

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

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

snuverink_j's avatar
snuverink_j committed
561
    double BcParameter[8] = {};
562

563 564
    BcParameter[0] = 0.001 * cycl_m->getRmin();
    BcParameter[1] = 0.001 * cycl_m->getRmax();
gsell's avatar
gsell committed
565

566
    // Store inner radius and outer radius of cyclotron field map in the list
567
    buildupFieldList(BcParameter, ElementBase::CYCLOTRON, cycl_m);
gsell's avatar
gsell committed
568 569 570 571 572 573 574 575 576 577
}

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

578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610

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
611 612 613 614 615
/**
 *
 *
 * @param coll
 */
616
void ParallelCyclotronTracker::visitCCollimator(const CCollimator &coll) {
gsell's avatar
gsell committed
617

618
    *gmsg << "* ------------------------------ Collimator ------------------------------" << endl;
gsell's avatar
gsell committed
619

620
    CCollimator* elptr = dynamic_cast<CCollimator *>(coll.clone());
621
    myElements.push_back(elptr);
gsell's avatar
gsell committed
622

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

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

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

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

635
    double zstart = elptr->getZStart();
636
    *gmsg << "* Zstart  = " << zstart << " [mm]" << endl;
637 638

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

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

frey_m's avatar
frey_m committed
644
    elptr->initialise(itsBunch_m);
gsell's avatar
gsell committed
645

snuverink_j's avatar
snuverink_j committed
646
    double BcParameter[8] = {};
647

648 649 650 651 652 653
    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 ;

654
    buildupFieldList(BcParameter, ElementBase::CCOLLIMATOR, elptr);
gsell's avatar
gsell committed
655 656 657 658 659 660 661 662 663 664 665 666
}

/**
 *
 *
 * @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
667 668 669 670 671 672 673 674 675 676 677 678
/**
 *
 *
 * @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
679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698
/**
 *
 *
 * @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()));
}

699 700
/**
 *
701 702
 *
 *  @param
703 704 705 706 707
 */
void ParallelCyclotronTracker::visitFlexibleCollimator(const FlexibleCollimator &) {

}

gsell's avatar
gsell committed
708 709 710 711 712 713 714 715 716 717
/**
 *
 *
 * @param lamb
 */
void ParallelCyclotronTracker::visitLambertson(const Lambertson &lamb) {
    *gmsg << "In Lambertson; L= " << lamb.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Lambertson *>(lamb.clone()));
}

718 719 720
void ParallelCyclotronTracker::visitOffset(const Offset & off) {
    if (opalRing_m == NULL)
        throw OpalException(
721 722
                            "ParallelCylcotronTracker::visitOffset",
                            "Attempt to place an offset when Ring not defined");
723 724
    Offset* offNonConst = const_cast<Offset*>(&off);
    offNonConst->updateGeometry(opalRing_m->getNextPosition(),
725
                                opalRing_m->getNextNormal());
726 727 728
    opalRing_m->appendElement(off);
}

gsell's avatar
gsell committed
729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760
/**
 *
 *
 * @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
761 762 763 764 765 766 767
/**
 *
 *
 * @param multT
 */
void ParallelCyclotronTracker::visitMultipoleT(const MultipoleT &multT) {
    *gmsg << "Adding MultipoleT" << endl;
768
    if (opalRing_m != NULL) {
ext-rogers_c's avatar
ext-rogers_c committed
769
        opalRing_m->appendElement(multT);
770
    } else {
ext-rogers_c's avatar
ext-rogers_c committed
771 772
        throw OpalException("ParallelCyclotronTracker::visitMultipoleT",
                            "Need to define a RINGDEFINITION to use MultipoleT element");
773
    }
ext-rogers_c's avatar
ext-rogers_c committed
774 775 776
    myElements.push_back(dynamic_cast<MultipoleT *>(multT.clone()));
}

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 816 817 818 819 820 821 822 823 824
/**
 *
 *
 * @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
825 826 827 828 829 830
/**
 *
 *
 * @param prob
 */
void ParallelCyclotronTracker::visitProbe(const Probe &prob) {
831
    *gmsg << "* ------------------------------  Probe ------------------------------" << endl;
832 833
    Probe *elptr = dynamic_cast<Probe *>(prob.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
834

835 836
    *gmsg << "* Name    = " << elptr->getName() << endl;

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

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

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

846
    double yend = elptr->getYEnd();
847
    *gmsg << "* YEnd    = " << yend << " [mm]" << endl;
gsell's avatar
gsell committed
848 849

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

snuverink_j's avatar
snuverink_j committed
852
    double BcParameter[8] = {};
853

854 855 856 857
    BcParameter[0] = 0.001 * xstart ;
    BcParameter[1] = 0.001 * xend;
    BcParameter[2] = 0.001 * ystart ;
    BcParameter[3] = 0.001 * yend;
858
    BcParameter[4] = 0.001 ; // width
gsell's avatar
gsell committed
859 860

    // store probe parameters in the list
861
    buildupFieldList(BcParameter, ElementBase::PROBE, elptr);
gsell's avatar
gsell committed
862 863 864 865 866 867 868 869 870 871 872 873
}

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

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

ext-rogers_c's avatar
ext-rogers_c committed
883 884
void ParallelCyclotronTracker::visitScalingFFAMagnet(const ScalingFFAMagnet &bend) {
    *gmsg << "Adding ScalingFFAMagnet" << endl;
ext-rogers_c's avatar
ext-rogers_c committed
885
    if (opalRing_m != NULL) {
886
        opalRing_m->appendElement(bend);
ext-rogers_c's avatar
ext-rogers_c committed
887
    } else {
ext-rogers_c's avatar
ext-rogers_c committed
888 889
        throw OpalException("ParallelCyclotronTracker::visitScalingFFAMagnet",
                            "Need to define a RINGDEFINITION to use ScalingFFAMagnet element");
ext-rogers_c's avatar
ext-rogers_c committed
890
    }
891 892
}

893 894 895 896 897 898
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",
899
                            "Need to define a RINGDEFINITION to use VariableRFCavity element");
900 901 902 903 904 905 906 907 908 909
}

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