ParallelCyclotronTracker.cpp 127 KB
Newer Older
gsell's avatar
gsell committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14
// ------------------------------------------------------------------------
// $RCSfile: ParallelCyclotronTracker.cpp,v $
// ------------------------------------------------------------------------
// $Revision: 1.1 $initialLocalNum_m
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: ParallelCyclotronTracker
//   The class for tracking particles with 3D space charge in Cyclotrons and FFAG's
//
// ------------------------------------------------------------------------
//
// $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/Ctunes.hpp>
kraus's avatar
kraus committed
20
#include "Algorithms/ParallelCyclotronTracker.h"
21 22
#include "Algorithms/PolynomialTimeDependence.h"
#include "Elements/OpalPolynomialTimeDependence.h"
23

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

63 64
#include "AbstractObjects/Element.h"

65
#include "Beamlines/FlaggedBeamline.h"
66
#include "Elements/OpalBeamline.h"
kraus's avatar
kraus committed
67
#include "AbsBeamline/Ring.h"
gsell's avatar
gsell committed
68 69 70

#include "BeamlineGeometry/Euclid3D.h"
#include "BeamlineGeometry/PlanarArcGeometry.h"
Jianjun Yang's avatar
Jianjun Yang committed
71
#include "BeamlineGeometry/RBendGeometry.h"
ext-rogers_c's avatar
ext-rogers_c committed
72
#include "BeamlineGeometry/StraightGeometry.h"
gsell's avatar
gsell committed
73 74 75 76 77 78 79 80 81 82
#include "Beamlines/Beamline.h"

#include "Fields/BMultipoleField.h"
#include "FixedAlgebra/FTps.h"
#include "FixedAlgebra/FTpsMath.h"
#include "FixedAlgebra/FVps.h"

#include "Physics/Physics.h"

#include "Utilities/OpalException.h"
83
#include "Utilities/Options.h"
gsell's avatar
gsell committed
84

85
#include "BasicActions/DumpFields.h"
86
#include "BasicActions/DumpEMFields.h"
87

88
#include "Structure/H5PartWrapperForPC.h"
89
#include "Structure/BoundaryGeometry.h"
gsell's avatar
gsell committed
90 91 92 93 94 95 96

#include "Ctunes.h"
#include <cassert>

#include <hdf5.h>
#include "H5hut.h"

frey_m's avatar
frey_m committed
97 98
//FIXME Remove headers and dynamic_cast in readOneBunchFromFile
#include "Algorithms/PartBunch.h"
99
#ifdef ENABLE_AMR
frey_m's avatar
frey_m committed
100 101
    #include "Algorithms/AmrPartBunch.h"
#endif
frey_m's avatar
frey_m committed
102

gsell's avatar
gsell committed
103 104
class Beamline;
class PartData;
105

gsell's avatar
gsell committed
106 107 108
using Physics::pi;
using Physics::q_e;

109

110
const double c_mmtns = Physics::c * 1.0e-6; // m/s --> mm/ns
111

112 113 114
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
115

snuverink_j's avatar
snuverink_j committed
116
//#define PSdim 6
gsell's avatar
gsell committed
117 118 119 120 121 122 123 124 125 126 127 128 129 130

extern Inform *gmsg;

// typedef FVector<double, PSdim> Vector;

/**
 * Constructor ParallelCyclotronTracker
 *
 * @param beamline
 * @param reference
 * @param revBeam
 * @param revTrack
 */
ParallelCyclotronTracker::ParallelCyclotronTracker(const Beamline &beamline,
131
                                                   const PartData &reference,
frey_m's avatar
frey_m committed
132 133 134 135 136
                                                   bool revBeam, bool revTrack)
    : Tracker(beamline, reference, revBeam, revTrack)
    , itsDataSink(nullptr)
    , itsMBDump_m(new MultiBunchDump())
    , bgf_m(nullptr)
137
    , numBunch_m(1)
frey_m's avatar
frey_m committed
138 139 140 141 142 143 144 145 146 147 148
    , lastDumpedStep_m(0)
    , eta_m(0.01)
    , myNode_m(Ippl::myNode())
    , initialLocalNum_m(0)
    , initialTotalNum_m(0)
    , onebunch_m(OpalData::getInstance()->getInputBasename() + "-onebunch.h5")
    , opalRing_m(nullptr)
    , itsStepper_mp(nullptr)
    , mode_m(MODE::UNDEFINED)
    , stepper_m(stepper::INTEGRATOR::UNDEFINED)
{
gsell's avatar
gsell committed
149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164
    itsBeamline = dynamic_cast<Beamline *>(beamline.clone());
}

/**
 * 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
165
                                                   PartBunchBase<double, 3> *bunch,
gsell's avatar
gsell committed
166 167 168
                                                   DataSink &ds,
                                                   const PartData &reference,
                                                   bool revBeam, bool revTrack,
frey_m's avatar
frey_m committed
169
                                                   int maxSTEPS, int timeIntegrator,
frey_m's avatar
frey_m committed
170
                                                   int numBunch)
frey_m's avatar
frey_m committed
171 172 173 174
    : Tracker(beamline, bunch, reference, revBeam, revTrack)
    , itsMBDump_m(new MultiBunchDump())
    , bgf_m(nullptr)
    , maxSteps_m(maxSTEPS)
175
    , numBunch_m(numBunch)
frey_m's avatar
frey_m committed
176 177 178 179 180 181 182 183 184
    , lastDumpedStep_m(0)
    , eta_m(0.01)
    , myNode_m(Ippl::myNode())
    , initialLocalNum_m(bunch->getLocalNum())
    , initialTotalNum_m(bunch->getTotalNum())
    , onebunch_m(OpalData::getInstance()->getInputBasename() + "-onebunch.h5")
    , opalRing_m(nullptr)
    , itsStepper_mp(nullptr)
{
gsell's avatar
gsell committed
185 186
    itsBeamline = dynamic_cast<Beamline *>(beamline.clone());
    itsDataSink = &ds;
frey_m's avatar
frey_m committed
187
    //  scaleFactor_m = itsBunch_m->getdT() * c;
gsell's avatar
gsell committed
188
    scaleFactor_m = 1;
189
    multiBunchMode_m = MB_MODE::NONE;
gsell's avatar
gsell committed
190 191 192 193 194

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

196 197 198 199 200 201 202 203 204
    // 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;
205

206 207 208 209 210 211 212 213
    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
214 215 216 217 218 219 220
}

/**
 * Destructor ParallelCyclotronTracker
 *
 */
ParallelCyclotronTracker::~ParallelCyclotronTracker() {
221 222
    for(Component* component : myElements) {
        delete(component);
gsell's avatar
gsell committed
223
    }
224 225
    for(auto fd : FieldDimensions) {
        delete(fd);
gsell's avatar
gsell committed
226 227
    }
    delete itsBeamline;
228
    // delete opalRing_m;
gsell's avatar
gsell committed
229 230
}

231 232 233 234 235 236
/**
 * AAA
 *
 * @param none
 */
void ParallelCyclotronTracker::initializeBoundaryGeometry() {
237 238
    for(Component * component : myElements) {
        bgf_m = dynamic_cast<ElementBase *>(component)->getBoundaryGeometry();
239 240 241 242 243 244 245 246 247 248
        if(!bgf_m)
            continue;
        else
            break;
    }
    if (bgf_m) {
        itsDataSink->writeGeomToVtk(*bgf_m, std::string("data/testGeometry-00000.vtk"));
        OpalData::getInstance()->setGlobalGeometry(bgf_m);
        *gmsg << "* Boundary geometry initialized " << endl;
    }
249
}
250

251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272

/// set the working sub-mode for multi-bunch mode: "FORCE" or "AUTO"
inline
void ParallelCyclotronTracker::setMultiBunchMode(const std::string& mbmode)
{
    if ( mbmode.compare("FORCE") == 0 ) {
        *gmsg << "FORCE mode: The multi bunches will be injected consecutively "
              << "after each revolution, until get \"TURNS\" bunches." << endl;
        multiBunchMode_m = MB_MODE::FORCE;
    } else if ( mbmode.compare("AUTO") == 0 ) {
        *gmsg << "AUTO mode: The multi bunches will be injected only when "
                      << "the distance between two neighboring bunches " << endl
                      << "is below the limitation. The control parameter is set to "
                      << CoeffDBunches_m << endl;
        multiBunchMode_m = MB_MODE::AUTO;
    } else if ( mbmode.compare("NONE") == 0 )
        multiBunchMode_m = MB_MODE::NONE;
    else
        throw OpalException("ParallelCyclotronTracker::setMultiBunchMode()",
                            "MBMODE name \"" + mbmode + "\" unknown.");
}

frey_m's avatar
frey_m committed
273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288
void ParallelCyclotronTracker::setMultiBunchBinning(std::string binning) {
    
    binning = Util::toUpper(binning);
    
    if ( binning.compare("BUNCH") == 0 ) {
        *gmsg << "Use 'BUNCH' injection for binnning." << endl;
        binningType_m = MB_BINNING::BUNCH;
    } else if ( binning.compare("GAMMA") == 0 ) {
        *gmsg << "Use 'GAMMA' for binning." << endl;
        binningType_m = MB_BINNING::GAMMA;
    } else {
        throw OpalException("ParallelCyclotronTracker::setMultiBunchBinning()",
                            "MB_BINNING name \"" + binning + "\" unknown.");
    }
}

289 290 291 292 293 294
/**
 *
 *
 * @param fn Base file name
 */
void ParallelCyclotronTracker::bgf_main_collision_test() {
295
    if(!bgf_m) return;
296

297
    Inform msg("bgf_main_collision_test ");
298

299 300 301
    /**
     *Here we check if a particles is outside the domain, flag it for deletion
     */
302

303
    Vector_t intecoords = 0.0;
304

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

308
    int triId = 0;
frey_m's avatar
frey_m committed
309 310 311
    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);
312
        if(res >= 0) {
frey_m's avatar
frey_m committed
313
            itsBunch_m->Bin[i] = -1;
314
        }
315
    }
316 317 318
}


gsell's avatar
gsell committed
319 320 321 322 323
/**
 *
 *
 * @param fn Base file name
 */
324
void ParallelCyclotronTracker::openFiles(std::string SfileName) {
gsell's avatar
gsell committed
325

326
    std::string  SfileName2 = SfileName + std::string("-Angle0.dat");
gsell's avatar
gsell committed
327
    outfTheta0_m.precision(8);
328
    outfTheta0_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
329
    outfTheta0_m.open(SfileName2.c_str());
330 331 332
    outfTheta0_m << "#  r [mm]      beta_r*gamma       "
                 << "theta [mm]      beta_theta*gamma        "
                 << "z [mm]          beta_z*gamma" << std::endl;
gsell's avatar
gsell committed
333

334
    SfileName2 = SfileName + std::string("-Angle1.dat");
gsell's avatar
gsell committed
335
    outfTheta1_m.precision(8);
336
    outfTheta1_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
337
    outfTheta1_m.open(SfileName2.c_str());
338 339 340
    outfTheta1_m << "#  r [mm]      beta_r*gamma       "
                 << "theta [mm]      beta_theta*gamma        "
                 << "z [mm]          beta_z*gamma"  << std::endl;
gsell's avatar
gsell committed
341

342
    SfileName2 = SfileName + std::string("-Angle2.dat");
gsell's avatar
gsell committed
343
    outfTheta2_m.precision(8);
344
    outfTheta2_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
345
    outfTheta2_m.open(SfileName2.c_str());
346 347 348
    outfTheta2_m << "#  r [mm]      beta_r*gamma       "
                 << "theta [mm]      beta_theta*gamma        "
                 << "z [mm]          beta_z*gamma"  << std::endl;
gsell's avatar
gsell committed
349 350 351

    // for single Particle Mode, output after each turn, to define matched initial phase ellipse.

352
    SfileName2 = SfileName + std::string("-afterEachTurn.dat");
gsell's avatar
gsell committed
353 354

    outfThetaEachTurn_m.precision(8);
355
    outfThetaEachTurn_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
356 357

    outfThetaEachTurn_m.open(SfileName2.c_str());
358 359 360
    outfTheta2_m << "#  r [mm]      beta_r*gamma       "
                 << "theta [mm]      beta_theta*gamma        "
                 << "z [mm]          beta_z*gamma"  << std::endl;
gsell's avatar
gsell committed
361 362 363 364 365 366 367 368 369 370 371 372 373 374 375
}

/**
 * Close all files related to
 * special output in the Cyclotron
 * mode.
 */
void ParallelCyclotronTracker::closeFiles() {

    outfTheta0_m.close();
    outfTheta1_m.close();
    outfTheta2_m.close();
    outfThetaEachTurn_m.close();
}

376
/**
377 378 379
 *
 * @param ring
 */
kraus's avatar
kraus committed
380
void ParallelCyclotronTracker::visitRing(const Ring &ring) {
Daniel Winklehner's avatar
Daniel Winklehner committed
381

382
    *gmsg << "* ----------------------------- Adding Ring ------------------------------ *" << endl;
Daniel Winklehner's avatar
Daniel Winklehner committed
383

snuverink_j's avatar
snuverink_j committed
384
    delete opalRing_m;
Daniel Winklehner's avatar
Daniel Winklehner committed
385

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

388
    myElements.push_back(opalRing_m);
Daniel Winklehner's avatar
Daniel Winklehner committed
389

frey_m's avatar
frey_m committed
390
    opalRing_m->initialise(itsBunch_m);
391 392 393 394

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

396
    if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
kraus's avatar
kraus committed
397
        throw OpalException("Error in ParallelCyclotronTracker::visitRing",
398 399
                            "PHIINIT is out of [-180, 180)!");
    }
Daniel Winklehner's avatar
Daniel Winklehner committed
400 401

    referenceZ = 0.0;
402
    referencePz = 0.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
403 404 405 406

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

407 408
    if(referencePtot < 0.0)
        referencePt *= -1.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
409

410 411
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
gsell's avatar
gsell committed
412

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

415
    buildupFieldList(BcParameter, ElementBase::RING, opalRing_m);
416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431

    // 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
432 433 434 435 436 437 438 439

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

440
    *gmsg << "* -------------------------- Adding Cyclotron ---------------------------- *" << endl;
gsell's avatar
gsell committed
441

442 443
    Cyclotron *elptr = dynamic_cast<Cyclotron *>(cycl.clone());
    myElements.push_back(elptr);
444

445
    // Is this a Spiral Inflector Simulation? If yes, we'll give the user some
446 447 448 449 450 451
    // useful information
    spiral_flag = elptr->getSpiralFlag();

    if(spiral_flag) {

        *gmsg << endl << "* This is a Spiral Inflector Simulation! This means the following:" << endl;
452
        *gmsg         << "* 1.) It is up to the user to provide appropriate geometry, electric and magnetic fields!" << endl;
453 454
        *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
455
        *gmsg         << "* 2.) For high currents it is strongly recommended to use the SAAMG fieldsolver," << endl;
456 457
        *gmsg         << "*     FFT does not give the correct results (boundaty conditions are missing)." << endl;
        *gmsg         << "* 3.) The whole geometry will be meshed and used for the fieldsolve." << endl;
458 459
        *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;
460 461 462 463 464 465 466
        *gmsg         << "*     (This is not an issue for spiral inflectors as they are typically < 100 keV/amu.)" << endl;
        *gmsg << endl << "* Note: For now, multi-bunch mode (MBM) needs to be de-activated for spiral inflector" << endl;
        *gmsg         << "* and space charge needs to be solved every time-step. numBunch_m and scSolveFreq are reset." << endl;
        numBunch_m = 1;

    }

467
    // Fresh run (no restart):
468
    if(!OpalData::getInstance()->inRestartRun()) {
469

470
        // Get reference values from cyclotron element
471
        // For now, these are still stored in mm. should be the only ones. -DW
472 473 474 475
        referenceR     = elptr->getRinit();
        referenceTheta = elptr->getPHIinit();
        referenceZ     = elptr->getZinit();
        referencePr    = elptr->getPRinit();
476 477
        referencePz    = elptr->getPZinit();

478
        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
479 480
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                "PHIINIT is out of [-180, 180)!");
481 482 483 484 485
        }

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

        // Calculate reference azimuthal (tangential) momentum from total-, z- and radial momentum:
486
        float insqrt = referencePtot * referencePtot - \
487
            referencePr * referencePr - referencePz * referencePz;
488 489 490 491 492

        if(insqrt < 0) {

            if(insqrt > -1.0e-10) {

493
                referencePt = 0.0;
494 495 496

            } else {

497
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
498
                                    "Pt imaginary!");
499 500 501 502 503 504 505
            }

        } else {

            referencePt = sqrt(insqrt);
        }

506
        if(referencePtot < 0.0)
Daniel Winklehner's avatar
Daniel Winklehner committed
507
            referencePt *= -1.0;
508 509
        // End calculate referencePt

510
        // Restart a run:
511 512
    } else {

513
        // If the user wants to save the restarted run in local frame,
514
        // make sure the previous h5 file was local too
515
      if (Options::psDumpFrame != Options::GLOBAL) {
winklehner_d's avatar
-DW  
winklehner_d committed
516
        if (!previousH5Local) {
517 518
                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
519
        }
520 521
            // 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
522 523
      } else {
        if (previousH5Local) {
524 525
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                    "You are trying a global restart from a local h5 file!");
526
            }
527
        }
528

529
        // Adjust some of the reference variables from the h5 file
530 531
        referencePhi *= Physics::deg2rad;
        referencePsi *= Physics::deg2rad;
532
        referencePtot = bega;
533
        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
534
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
535 536
                                "PHIINIT is out of [-180, 180)!");
        }
537 538
    }

539 540
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
541

542
    *gmsg << endl;
adelmann's avatar
adelmann committed
543
    *gmsg << "* Bunch global starting position:" << endl;
544 545
    *gmsg << "* RINIT = " << referenceR  << " [mm]" << endl;
    *gmsg << "* PHIINIT = " << referenceTheta << " [deg]" << endl;
546
    *gmsg << "* ZINIT = " << referenceZ << " [mm]" << endl;
547
    *gmsg << endl;
adelmann's avatar
adelmann committed
548
    *gmsg << "* Bunch global starting momenta:" << endl;
gsell's avatar
gsell committed
549 550
    *gmsg << "* Initial gamma = " << itsReference.getGamma() << endl;
    *gmsg << "* Initial beta = " << itsReference.getBeta() << endl;
551
    *gmsg << "* Reference total momentum (beta * gamma) = " << referencePtot * 1000.0 << " [MCU]" << endl;
552 553 554
    *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;
555
    *gmsg << endl;
adelmann's avatar
adelmann committed
556

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

560 561 562
    // 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
563

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

567
    std::string type = elptr->getCyclotronType();
568
    *gmsg << "* Type of cyclotron = " << type << " " << endl;
569

570 571
    double rmin = elptr->getMinR();
    double rmax = elptr->getMaxR();
adelmann's avatar
adelmann committed
572
    *gmsg << "* Radial aperture = " << rmin << " ... " << rmax<<" [m] "<< endl;
573 574 575

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

    double h = elptr->getCyclHarm();
579
    *gmsg << "* Number of trimcoils = " << elptr->getNumberOfTrimcoils() << endl;
580
    *gmsg << "* Harmonic number h = " << h << " " << endl;
581

gsell's avatar
gsell committed
582 583 584 585 586 587 588 589 590
    /**
     * 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
     * fieldflag = 5, readin FFAG format file for MSU/FNAL FFAG
     * fieldflag = 6, readin both median plane B field map and 3D E field map of RF cavity for compact cyclotron
591
     * fieldflag = 7, read in fields for Daniel's synchrocyclotron simulations
gsell's avatar
gsell committed
592
     */
593
    int fieldflag = elptr->getFieldFlag(type);
gsell's avatar
gsell committed
594

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

snuverink_j's avatar
snuverink_j committed
598
    double BcParameter[8] = {};
599

600 601
    BcParameter[0] = 0.001 * elptr->getRmin();
    BcParameter[1] = 0.001 * elptr->getRmax();
gsell's avatar
gsell committed
602

603
    // Store inner radius and outer radius of cyclotron field map in the list
604
    buildupFieldList(BcParameter, ElementBase::CYCLOTRON, elptr);
gsell's avatar
gsell committed
605 606 607 608 609 610 611 612 613 614 615 616 617 618 619
}

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

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

622
    *gmsg << "* --------- Collimator -----------------------------" << endl;
gsell's avatar
gsell committed
623

624
    CCollimator* elptr = dynamic_cast<CCollimator *>(coll.clone());
625
    myElements.push_back(elptr);
gsell's avatar
gsell committed
626

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

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

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

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

639
    double zstart = elptr->getZStart();
adelmann's avatar
adelmann committed
640
    *gmsg << "* Zstart= " << zstart << " [mm]" << endl;
641 642

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

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

frey_m's avatar
frey_m committed
648
    elptr->initialise(itsBunch_m);
gsell's avatar
gsell committed
649

snuverink_j's avatar
snuverink_j committed
650
    double BcParameter[8] = {};
651

652 653 654 655 656 657
    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 ;

658
    buildupFieldList(BcParameter, ElementBase::CCOLLIMATOR, elptr);
gsell's avatar
gsell committed
659 660 661 662 663 664 665 666 667 668 669 670
}

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

703 704
/**
 *
705 706
 *
 *  @param
707 708 709 710 711
 */
void ParallelCyclotronTracker::visitFlexibleCollimator(const FlexibleCollimator &) {

}

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

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

gsell's avatar
gsell committed
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
/**
 *
 *
 * @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
765 766 767 768 769 770 771
/**
 *
 *
 * @param multT
 */
void ParallelCyclotronTracker::visitMultipoleT(const MultipoleT &multT) {
    *gmsg << "Adding MultipoleT" << endl;
772
    if (opalRing_m != NULL) {
ext-rogers_c's avatar
ext-rogers_c committed
773
        opalRing_m->appendElement(multT);
774
    } else {
ext-rogers_c's avatar
ext-rogers_c committed
775 776
        throw OpalException("ParallelCyclotronTracker::visitMultipoleT",
                            "Need to define a RINGDEFINITION to use MultipoleT element");
777
    }
ext-rogers_c's avatar
ext-rogers_c committed
778 779 780
    myElements.push_back(dynamic_cast<MultipoleT *>(multT.clone()));
}

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 825 826 827 828
/**
 *
 *
 * @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
829 830 831 832 833 834
/**
 *
 *
 * @param prob
 */
void ParallelCyclotronTracker::visitProbe(const Probe &prob) {
835
    *gmsg << "* -----------  Probe -------------------------------" << endl;
836 837
    Probe *elptr = dynamic_cast<Probe *>(prob.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
838

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

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

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

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

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

snuverink_j's avatar
snuverink_j committed
854
    double BcParameter[8] = {};
855

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

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

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

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

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

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

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

gsell's avatar
gsell committed
914 915 916 917 918 919 920 921
/**
 *
 *
 * @param as
 */
void ParallelCyclotronTracker::visitRFCavity(const RFCavity &as) {

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

923 924
    RFCavity *elptr = dynamic_cast<RFCavity *>(as.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
925 926 927 928 929 930 931

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

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

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

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

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

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

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

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

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

956 957 958
    /*
      Setup time dependence and in case of no
      timedependence use a polynom with  a_0 = 1 and a_k = 0, k = 1,2,3.
959
    */
960

961 962 963 964
    std::shared_ptr<AbstractTimeDependence> freq_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> ampl_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> phase_atd = nullptr;

965
    dvector_t  unityVec;
966 967 968 969
    unityVec.push_back(1.);