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

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

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

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

#include "Physics/Physics.h"

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

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

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

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

87

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

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

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

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

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

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

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

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

frey_m's avatar
frey_m committed
183

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

187
    Inform msg("bgf_main_collision_test ");
188

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

193
    Vector_t intecoords = 0.0;
194

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

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

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


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

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

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

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

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

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

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


283

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    if(spiral_flag) {

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

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

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

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

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

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

        if(insqrt < 0) {

            if(insqrt > -1.0e-10) {

439
                referencePt = 0.0;
440 441 442

            } else {

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

        } else {

            referencePt = sqrt(insqrt);
        }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

568
    *gmsg << "* --------- Collimator -----------------------------" << endl;
gsell's avatar
gsell committed
569

570
    CCollimator* elptr = dynamic_cast<CCollimator *>(coll.clone());
571
    myElements.push_back(elptr);
gsell's avatar
gsell committed
572

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

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

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

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

585
    double zstart = elptr->getZStart();
adelmann's avatar
adelmann committed
586
    *gmsg << "* Zstart= " << zstart << " [mm]" << endl;
587 588

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

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

frey_m's avatar
frey_m committed
594
    elptr->initialise(itsBunch_m);
gsell's avatar
gsell committed
595

snuverink_j's avatar
snuverink_j committed
596
    double BcParameter[8] = {};
597

598 599 600 601 602 603
    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 ;

604
    buildupFieldList(BcParameter, ElementBase::CCOLLIMATOR, elptr);
gsell's avatar
gsell committed
605 606 607 608 609 610 611 612 613 614 615 616
}

/**
 *
 *
 * @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
617 618 619 620 621 622 623 624 625 626 627 628
/**
 *
 *
 * @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
629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648
/**
 *
 *
 * @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()));
}

649 650
/**
 *
651 652
 *
 *  @param
653 654 655 656 657
 */
void ParallelCyclotronTracker::visitFlexibleCollimator(const FlexibleCollimator &) {

}

gsell's avatar
gsell committed
658 659 660 661 662 663 664 665 666 667
/**
 *
 *
 * @param lamb
 */
void ParallelCyclotronTracker::visitLambertson(const Lambertson &lamb) {
    *gmsg << "In Lambertson; L= " << lamb.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Lambertson *>(lamb.clone()));
}

668 669 670
void ParallelCyclotronTracker::visitOffset(const Offset & off) {
    if (opalRing_m == NULL)
        throw OpalException(
671 672
                            "ParallelCylcotronTracker::visitOffset",
                            "Attempt to place an offset when Ring not defined");
673 674
    Offset* offNonConst = const_cast<Offset*>(&off);
    offNonConst->updateGeometry(opalRing_m->getNextPosition(),
675
                                opalRing_m->getNextNormal());
676 677 678
    opalRing_m->appendElement(off);
}

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 699 700 701 702 703 704 705 706 707 708 709 710
/**
 *
 *
 * @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
711 712 713 714 715 716 717
/**
 *
 *
 * @param multT
 */
void ParallelCyclotronTracker::visitMultipoleT(const MultipoleT &multT) {
    *gmsg << "Adding MultipoleT" << endl;
718
    if (opalRing_m != NULL) {
ext-rogers_c's avatar
ext-rogers_c committed
719
        opalRing_m->appendElement(multT);
720
    } else {
ext-rogers_c's avatar
ext-rogers_c committed
721 722
        throw OpalException("ParallelCyclotronTracker::visitMultipoleT",
                            "Need to define a RINGDEFINITION to use MultipoleT element");
723
    }
ext-rogers_c's avatar
ext-rogers_c committed
724 725 726
    myElements.push_back(dynamic_cast<MultipoleT *>(multT.clone()));
}

727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774
/**
 *
 *
 * @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
775 776 777 778 779 780
/**
 *
 *
 * @param prob
 */
void ParallelCyclotronTracker::visitProbe(const Probe &prob) {
781
    *gmsg << "* -----------  Probe -------------------------------" << endl;
782 783
    Probe *elptr = dynamic_cast<Probe *>(prob.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
784

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

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

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

794
    double yend = elptr->getYEnd();
795
    *gmsg << "YEnd= " << yend << " [mm]" << endl;
gsell's avatar
gsell committed
796 797

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

snuverink_j's avatar
snuverink_j committed
800
    double BcParameter[8] = {};
801

802 803 804 805
    BcParameter[0] = 0.001 * xstart ;
    BcParameter[1] = 0.001 * xend;
    BcParameter[2] = 0.001 * ystart ;
    BcParameter[3] = 0.001 * yend;
806
    BcParameter[4] = 0.001 ; // width
gsell's avatar
gsell committed
807 808

    // store probe parameters in the list
809
    buildupFieldList(BcParameter, ElementBase::PROBE, elptr);
gsell's avatar
gsell committed
810 811 812 813 814 815 816 817 818 819 820 821
}

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

822 823 824 825 826 827
void ParallelCyclotronTracker::visitSBend3D(const SBend3D &bend) {
    *gmsg << "Adding SBend3D" << endl;
    if (opalRing_m != NULL)
        opalRing_m->appendElement(bend);
    else
        throw OpalException("ParallelCyclotronTracker::visitSBend3D",
828
                            "Need to define a RINGDEFINITION to use SBend3D element");
829 830
}

ext-rogers_c's avatar
ext-rogers_c committed
831 832
void ParallelCyclotronTracker::visitScalingFFAMagnet(const ScalingFFAMagnet &bend) {
    *gmsg << "Adding ScalingFFAMagnet" << endl;
ext-rogers_c's avatar
ext-rogers_c committed
833
    if (opalRing_m != NULL) {
834
        opalRing_m->appendElement(bend);
ext-rogers_c's avatar
ext-rogers_c committed
835
    } else {
ext-rogers_c's avatar
ext-rogers_c committed
836 837
        throw OpalException("ParallelCyclotronTracker::visitScalingFFAMagnet",
                            "Need to define a RINGDEFINITION to use ScalingFFAMagnet element");
ext-rogers_c's avatar
ext-rogers_c committed
838
    }
839 840
}

841 842 843 844 845 846
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",
847
                            "Need to define a RINGDEFINITION to use VariableRFCavity element");
848 849 850 851 852 853 854 855 856 857
}

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

gsell's avatar
gsell committed
860 861 862 863 864 865 866 867
/**
 *
 *
 * @param as
 */
void ParallelCyclotronTracker::visitRFCavity(const RFCavity &as) {

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

869 870
    RFCavity *elptr = dynamic_cast<RFCavity *>(as.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
871 872 873 874 875 876 877

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

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

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

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

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

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

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

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

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

902 903 904
    /*
      Setup time dependence and in case of no
      timedependence use a polynom with  a_0 = 1 and a_k = 0, k = 1,2,3.
Christof Metzger-Kraus's avatar