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 27

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

30
#include "AbsBeamline/BeamStripping.h"
31
#include "AbsBeamline/CCollimator.h"
gsell's avatar
gsell committed
32 33
#include "AbsBeamline/Corrector.h"
#include "AbsBeamline/Cyclotron.h"
34
#include "AbsBeamline/CyclotronValley.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

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

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

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

#include "Physics/Physics.h"

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

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

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

gsell's avatar
gsell committed
85 86 87
using Physics::pi;
using Physics::q_e;

88

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

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

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

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

gsell's avatar
gsell committed
140 141 142 143
    IntegrationTimer_m = IpplTimings::getTimer("Integration");
    TransformTimer_m   = IpplTimings::getTimer("Frametransform");
    DumpTimer_m        = IpplTimings::getTimer("Dump");
    BinRepartTimer_m   = IpplTimings::getTimer("Binaryrepart");
144 145
    PluginElemTimer_m  = IpplTimings::getTimer("PluginElements");
    DelParticleTimer_m = IpplTimings::getTimer("DeleteParticles");
146

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

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

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

frey_m's avatar
frey_m committed
184

185
void ParallelCyclotronTracker::bgf_main_collision_test() {
186
    if(!bgf_m) return;
187

188
    Inform msg("bgf_main_collision_test ");
189

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

194
    Vector_t intecoords = 0.0;
195

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

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

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


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

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

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

frey_m's avatar
frey_m committed
259 260 261
        // 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
262

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

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

frey_m's avatar
frey_m committed
276 277 278 279 280
    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
    }
281 282 283
}


284

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

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

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

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

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

328
    *gmsg << "* ----------------------------- Adding Ring ------------------------------ *" << endl;
Daniel Winklehner's avatar
Daniel Winklehner committed
329

snuverink_j's avatar
snuverink_j committed
330
    delete opalRing_m;
Daniel Winklehner's avatar
Daniel Winklehner committed
331

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

334
    myElements.push_back(opalRing_m);
Daniel Winklehner's avatar
Daniel Winklehner committed
335

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

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

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

    referenceZ = 0.0;
348
    referencePz = 0.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
349 350 351 352

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

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

356 357
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
gsell's avatar
gsell committed
358

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

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

    // 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
378 379 380 381 382 383 384 385

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

386
    *gmsg << "* -------------------------- Adding Cyclotron ---------------------------- *" << endl;
gsell's avatar
gsell committed
387

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

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

    if(spiral_flag) {

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

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

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

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

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

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

        if(insqrt < 0) {

            if(insqrt > -1.0e-10) {

440
                referencePt = 0.0;
441 442 443

            } else {

444
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
445
                                    "Pt imaginary!");
446 447 448 449 450 451 452
            }

        } else {

            referencePt = sqrt(insqrt);
        }

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

457
        // Restart a run:
458 459
    } else {

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

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

486 487
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
488

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

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

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

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

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

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

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

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

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

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

snuverink_j's avatar
snuverink_j committed
545
    double BcParameter[8] = {};
546

547 548
    BcParameter[0] = 0.001 * cycl_m->getRmin();
    BcParameter[1] = 0.001 * cycl_m->getRmax();
gsell's avatar
gsell committed
549

550
    // Store inner radius and outer radius of cyclotron field map in the list
551
    buildupFieldList(BcParameter, ElementBase::CYCLOTRON, cycl_m);
gsell's avatar
gsell committed
552 553 554 555 556 557 558 559 560 561
}

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

562 563 564 565 566 567 568 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

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
595 596 597 598 599
/**
 *
 *
 * @param coll
 */
600
void ParallelCyclotronTracker::visitCCollimator(const CCollimator &coll) {
gsell's avatar
gsell committed
601

602
    *gmsg << "* ------------------------------ Collimator ------------------------------" << endl;
gsell's avatar
gsell committed
603

604
    CCollimator* elptr = dynamic_cast<CCollimator *>(coll.clone());
605
    myElements.push_back(elptr);
gsell's avatar
gsell committed
606

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

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

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

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

619
    double zstart = elptr->getZStart();
620
    *gmsg << "* Zstart  = " << zstart << " [mm]" << endl;
621 622

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

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

frey_m's avatar
frey_m committed
628
    elptr->initialise(itsBunch_m);
gsell's avatar
gsell committed
629

snuverink_j's avatar
snuverink_j committed
630
    double BcParameter[8] = {};
631

632 633 634 635 636 637
    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 ;

638
    buildupFieldList(BcParameter, ElementBase::CCOLLIMATOR, elptr);
gsell's avatar
gsell committed
639 640 641 642 643 644 645 646 647 648 649 650
}

/**
 *
 *
 * @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
651 652 653 654 655 656 657 658 659 660 661 662
/**
 *
 *
 * @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
663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682
/**
 *
 *
 * @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()));
}

683 684
/**
 *
685 686
 *
 *  @param
687 688 689 690 691
 */
void ParallelCyclotronTracker::visitFlexibleCollimator(const FlexibleCollimator &) {

}

gsell's avatar
gsell committed
692 693 694 695 696 697 698 699 700 701
/**
 *
 *
 * @param lamb
 */
void ParallelCyclotronTracker::visitLambertson(const Lambertson &lamb) {
    *gmsg << "In Lambertson; L= " << lamb.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Lambertson *>(lamb.clone()));
}

702 703 704
void ParallelCyclotronTracker::visitOffset(const Offset & off) {
    if (opalRing_m == NULL)
        throw OpalException(
705 706
                            "ParallelCylcotronTracker::visitOffset",
                            "Attempt to place an offset when Ring not defined");
707 708
    Offset* offNonConst = const_cast<Offset*>(&off);
    offNonConst->updateGeometry(opalRing_m->getNextPosition(),
709
                                opalRing_m->getNextNormal());
710 711 712
    opalRing_m->appendElement(off);
}

gsell's avatar
gsell committed
713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744
/**
 *
 *
 * @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
745 746 747 748 749 750 751
/**
 *
 *
 * @param multT
 */
void ParallelCyclotronTracker::visitMultipoleT(const MultipoleT &multT) {
    *gmsg << "Adding MultipoleT" << endl;
752
    if (opalRing_m != NULL) {
ext-rogers_c's avatar
ext-rogers_c committed
753
        opalRing_m->appendElement(multT);
754
    } else {
ext-rogers_c's avatar
ext-rogers_c committed
755 756
        throw OpalException("ParallelCyclotronTracker::visitMultipoleT",
                            "Need to define a RINGDEFINITION to use MultipoleT element");
757
    }
ext-rogers_c's avatar
ext-rogers_c committed
758 759 760
    myElements.push_back(dynamic_cast<MultipoleT *>(multT.clone()));
}

761 762 763 764 765 766 767 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
/**
 *
 *
 * @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
809 810 811 812 813 814
/**
 *
 *
 * @param prob
 */
void ParallelCyclotronTracker::visitProbe(const Probe &prob) {
815
    *gmsg << "* ------------------------------  Probe ------------------------------" << endl;
816 817
    Probe *elptr = dynamic_cast<Probe *>(prob.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
818

819 820
    *gmsg << "* Name    = " << elptr->getName() << endl;

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

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

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

830
    double yend = elptr->getYEnd();
831
    *gmsg << "* YEnd    = " << yend << " [mm]" << endl;
gsell's avatar
gsell committed
832 833

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

snuverink_j's avatar
snuverink_j committed
836
    double BcParameter[8] = {};
837

838 839 840 841
    BcParameter[0] = 0.001 * xstart ;
    BcParameter[1] = 0.001 * xend;
    BcParameter[2] = 0.001 * ystart ;
    BcParameter[3] = 0.001 * yend;
842
    BcParameter[4] = 0.001 ; // width
gsell's avatar
gsell committed
843 844

    // store probe parameters in the list
845
    buildupFieldList(BcParameter, ElementBase::PROBE, elptr);
gsell's avatar
gsell committed
846 847 848 849 850 851 852 853 854 855 856 857
}

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

858 859 860 861 862 863
void ParallelCyclotronTracker::visitSBend3D(const SBend3D &bend) {
    *gmsg << "Adding SBend3D" << endl;
    if (opalRing_m != NULL)
        opalRing_m->appendElement(bend);
    else
        throw OpalException("ParallelCyclotronTracker::visitSBend3D",
864
                            "Need to define a RINGDEFINITION to use SBend3D element");
865 866
}

ext-rogers_c's avatar
ext-rogers_c committed
867 868
void ParallelCyclotronTracker::visitScalingFFAMagnet(const ScalingFFAMagnet &bend) {
    *gmsg << "Adding ScalingFFAMagnet" << endl;
ext-rogers_c's avatar
ext-rogers_c committed
869
    if (opalRing_m != NULL) {
870
        opalRing_m->appendElement(bend);
ext-rogers_c's avatar
ext-rogers_c committed
871
    } else {
ext-rogers_c's avatar
ext-rogers_c committed
872 873
        throw OpalException("ParallelCyclotronTracker::visitScalingFFAMagnet",
                            "Need to define a RINGDEFINITION to use ScalingFFAMagnet element");
ext-rogers_c's avatar
ext-rogers_c committed
874
    }
875 876
}

877 878 879 880 881 882
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",
883
                            "Need to define a RINGDEFINITION to use VariableRFCavity element");
884 885 886 887 888 889 890 891 892 893
}

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

gsell's avatar
gsell committed
896 897 898 899 900 901 902
/**
 *
 *
 * @param as
 */
void ParallelCyclotronTracker::visitRFCavity(const RFCavity &as) {

903
    *gmsg << "* ------------------------------ RFCavity ------------------------------" << endl;
904

905 906
    RFCavity *elptr = dynamic_cast<RFCavity *>(as.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
907 908 909 910 911 912 913

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

914
    double rmin = elptr->getRmin();
gsell's avatar