ParallelCyclotronTracker.cpp 126 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 31 32

#include "AbsBeamline/Collimator.h"
#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 42 43 44 45 46
#include "AbsBeamline/Marker.h"
#include "AbsBeamline/Monitor.h"
#include "AbsBeamline/Multipole.h"
#include "AbsBeamline/Probe.h"
#include "AbsBeamline/RBend.h"
#include "AbsBeamline/RFCavity.h"
#include "AbsBeamline/RFQuadrupole.h"
#include "AbsBeamline/SBend.h"
47
#include "AbsBeamline/SBend3D.h"
48
#include "AbsBeamline/ScalingFFAGMagnet.h"
gsell's avatar
gsell committed
49 50 51 52 53
#include "AbsBeamline/Separator.h"
#include "AbsBeamline/Septum.h"
#include "AbsBeamline/Solenoid.h"
#include "AbsBeamline/CyclotronValley.h"
#include "AbsBeamline/Stripper.h"
54
#include "AbsBeamline/VariableRFCavity.h"
55

56 57
#include "AbstractObjects/Element.h"

58
#include "Elements/OpalBeamline.h"
kraus's avatar
kraus committed
59
#include "AbsBeamline/Ring.h"
gsell's avatar
gsell committed
60 61 62

#include "BeamlineGeometry/Euclid3D.h"
#include "BeamlineGeometry/PlanarArcGeometry.h"
Jianjun Yang's avatar
Jianjun Yang committed
63
#include "BeamlineGeometry/RBendGeometry.h"
gsell's avatar
gsell committed
64 65 66 67 68 69 70 71 72 73 74 75
#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/NumToStr.h"
#include "Utilities/OpalException.h"

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

79
#include "Structure/H5PartWrapperForPC.h"
80
#include "Structure/BoundaryGeometry.h"
81
#include "Utilities/Options.h"
gsell's avatar
gsell committed
82 83 84 85 86 87 88

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

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

frey_m's avatar
frey_m committed
89 90
//FIXME Remove headers and dynamic_cast in readOneBunchFromFile
#include "Algorithms/PartBunch.h"
91
#ifdef ENABLE_AMR
frey_m's avatar
frey_m committed
92 93
    #include "Algorithms/AmrPartBunch.h"
#endif
frey_m's avatar
frey_m committed
94

gsell's avatar
gsell committed
95 96
class Beamline;
class PartData;
97

gsell's avatar
gsell committed
98 99 100
using Physics::pi;
using Physics::q_e;

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

103 104 105
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
106

snuverink_j's avatar
snuverink_j committed
107
//#define PSdim 6
gsell's avatar
gsell committed
108 109 110 111 112 113 114 115 116 117 118 119 120 121

extern Inform *gmsg;

// typedef FVector<double, PSdim> Vector;

/**
 * Constructor ParallelCyclotronTracker
 *
 * @param beamline
 * @param reference
 * @param revBeam
 * @param revTrack
 */
ParallelCyclotronTracker::ParallelCyclotronTracker(const Beamline &beamline,
122 123
                                                   const PartData &reference,
                                                   bool revBeam, bool revTrack):
gsell's avatar
gsell committed
124
    Tracker(beamline, reference, revBeam, revTrack),
125 126
    itsDataSink(nullptr),
    bgf_m(nullptr),
127
    lastDumpedStep_m(0),
128
    eta_m(0.01),
129 130
    initialR_m(nullptr),
    initialP_m(nullptr),
gsell's avatar
gsell committed
131 132
    myNode_m(Ippl::myNode()),
    initialLocalNum_m(0),
133
    initialTotalNum_m(0),
134
    opalRing_m(nullptr),
135 136
    itsStepper_mp(nullptr),
    mode_m(MODE::UNDEFINED),
137
    stepper_m(stepper::INTEGRATOR::UNDEFINED) {
gsell's avatar
gsell committed
138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153
    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
154
                                                   PartBunchBase<double, 3> *bunch,
gsell's avatar
gsell committed
155 156 157
                                                   DataSink &ds,
                                                   const PartData &reference,
                                                   bool revBeam, bool revTrack,
158
                                                   int maxSTEPS, int timeIntegrator):
159
    Tracker(beamline, bunch, reference, revBeam, revTrack),
160
    bgf_m(nullptr),
gsell's avatar
gsell committed
161
    maxSteps_m(maxSTEPS),
162
    lastDumpedStep_m(0),
163
    eta_m(0.01),
164 165
    initialR_m(nullptr),
    initialP_m(nullptr),
gsell's avatar
gsell committed
166
    myNode_m(Ippl::myNode()),
frey_m's avatar
frey_m committed
167 168
    initialLocalNum_m(bunch->getLocalNum()),
    initialTotalNum_m(bunch->getTotalNum()),
169 170
    opalRing_m(nullptr),
    itsStepper_mp(nullptr) {
gsell's avatar
gsell committed
171 172
    itsBeamline = dynamic_cast<Beamline *>(beamline.clone());
    itsDataSink = &ds;
frey_m's avatar
frey_m committed
173
    //  scaleFactor_m = itsBunch_m->getdT() * c;
gsell's avatar
gsell committed
174 175 176 177 178 179 180
    scaleFactor_m = 1;
    multiBunchMode_m = 0;

    IntegrationTimer_m = IpplTimings::getTimer("Integration");
    TransformTimer_m   = IpplTimings::getTimer("Frametransform");
    DumpTimer_m        = IpplTimings::getTimer("Dump");
    BinRepartTimer_m   = IpplTimings::getTimer("Binaryrepart");
181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199
    
    // 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;
    
    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
200 201 202 203 204 205 206
}

/**
 * Destructor ParallelCyclotronTracker
 *
 */
ParallelCyclotronTracker::~ParallelCyclotronTracker() {
207 208
    for(Component* component : myElements) {
        delete(component);
gsell's avatar
gsell committed
209
    }
210 211
    for(auto fd : FieldDimensions) {
        delete(fd);
gsell's avatar
gsell committed
212 213
    }
    delete itsBeamline;
214
    // delete opalRing_m;
gsell's avatar
gsell committed
215 216
}

217 218 219 220 221 222
/**
 * AAA
 *
 * @param none
 */
void ParallelCyclotronTracker::initializeBoundaryGeometry() {
223 224
    for(Component * component : myElements) {
        bgf_m = dynamic_cast<ElementBase *>(component)->getBoundaryGeometry();
225 226 227 228 229 230 231 232 233 234
        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;
    }
235
}
236

237 238 239 240 241 242
/**
 *
 *
 * @param fn Base file name
 */
void ParallelCyclotronTracker::bgf_main_collision_test() {
243
    if(!bgf_m) return;
244

245
    Inform msg("bgf_main_collision_test ");
246

247 248 249
    /**
     *Here we check if a particles is outside the domain, flag it for deletion
     */
250

251
    Vector_t intecoords = 0.0;
252

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

256
    int triId = 0;
frey_m's avatar
frey_m committed
257 258 259
    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);
260
        if(res >= 0) {
frey_m's avatar
frey_m committed
261
            itsBunch_m->Bin[i] = -1;
262
        }
263
    }
264 265 266
}


gsell's avatar
gsell committed
267 268 269 270 271
/**
 *
 *
 * @param fn Base file name
 */
272
void ParallelCyclotronTracker::openFiles(std::string SfileName) {
gsell's avatar
gsell committed
273

274
    std::string  SfileName2 = SfileName + std::string("-Angle0.dat");
gsell's avatar
gsell committed
275
    outfTheta0_m.precision(8);
276
    outfTheta0_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
277
    outfTheta0_m.open(SfileName2.c_str());
278 279 280
    outfTheta0_m << "#  r [mm]      beta_r*gamma       "
                 << "theta [mm]      beta_theta*gamma        "
                 << "z [mm]          beta_z*gamma" << std::endl;
gsell's avatar
gsell committed
281

282
    SfileName2 = SfileName + std::string("-Angle1.dat");
gsell's avatar
gsell committed
283
    outfTheta1_m.precision(8);
284
    outfTheta1_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
285
    outfTheta1_m.open(SfileName2.c_str());
286 287 288
    outfTheta1_m << "#  r [mm]      beta_r*gamma       "
                 << "theta [mm]      beta_theta*gamma        "
                 << "z [mm]          beta_z*gamma"  << std::endl;
gsell's avatar
gsell committed
289

290
    SfileName2 = SfileName + std::string("-Angle2.dat");
gsell's avatar
gsell committed
291
    outfTheta2_m.precision(8);
292
    outfTheta2_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
293
    outfTheta2_m.open(SfileName2.c_str());
294 295 296
    outfTheta2_m << "#  r [mm]      beta_r*gamma       "
                 << "theta [mm]      beta_theta*gamma        "
                 << "z [mm]          beta_z*gamma"  << std::endl;
gsell's avatar
gsell committed
297 298 299

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

300
    SfileName2 = SfileName + std::string("-afterEachTurn.dat");
gsell's avatar
gsell committed
301 302

    outfThetaEachTurn_m.precision(8);
303
    outfThetaEachTurn_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
304 305

    outfThetaEachTurn_m.open(SfileName2.c_str());
306 307 308
    outfTheta2_m << "#  r [mm]      beta_r*gamma       "
                 << "theta [mm]      beta_theta*gamma        "
                 << "z [mm]          beta_z*gamma"  << std::endl;
gsell's avatar
gsell committed
309 310 311 312 313 314 315 316 317 318 319 320 321 322 323
}

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

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

330
    *gmsg << "* ----------------------------- Adding Ring ------------------------------ *" << endl;
Daniel Winklehner's avatar
Daniel Winklehner committed
331 332

    if (opalRing_m != NULL)
333
        delete opalRing_m;
Daniel Winklehner's avatar
Daniel Winklehner committed
334

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

337
    myElements.push_back(opalRing_m);
Daniel Winklehner's avatar
Daniel Winklehner committed
338

frey_m's avatar
frey_m committed
339
    opalRing_m->initialise(itsBunch_m);
340 341 342 343

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

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

    referenceZ = 0.0;
351
    referencePz = 0.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
352 353 354 355

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

356 357
    if(referencePtot < 0.0)
        referencePt *= -1.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
358

359 360
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
gsell's avatar
gsell committed
361

362
    double BcParameter[8];
Daniel Winklehner's avatar
Daniel Winklehner committed
363

364
    for(int i = 0; i < 8; i++)
Daniel Winklehner's avatar
Daniel Winklehner committed
365 366
        BcParameter[i] = 0.0;

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

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

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

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

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

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

    if(spiral_flag) {

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

    }

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

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

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

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

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

        if(insqrt < 0) {

            if(insqrt > -1.0e-10) {

	        referencePt = 0.0;

            } else {

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

        } else {

            referencePt = sqrt(insqrt);
        }

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

462
        // Restart a run:
463 464
    } else {

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

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

491 492
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
493

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

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

512 513 514
    // 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
515

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

519
    std::string type = elptr->getCyclotronType();
520
    *gmsg << "* Type of cyclotron = " << type << " " << endl;
521

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

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

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

gsell's avatar
gsell committed
534 535 536 537 538 539 540 541 542
    /**
     * To ease the initialise() function, set a integral parameter fieldflag internally.
     * Its value is  by the option "TYPE" of the element  "CYCLOTRON"
     * fieldflag = 1, readin PSI format measured field file (default)
     * fieldflag = 2, readin carbon cyclotron field file created by Jianjun Yang, TYPE=CARBONCYCL
     * fieldflag = 3, readin ANSYS format file for CYCIAE-100 created by Jianjun Yang, TYPE=CYCIAE
     * fieldflag = 4, readin AVFEQ format file for Riken cyclotrons
     * 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
543
     * fieldflag = 7, read in fields for Daniel's synchrocyclotron simulations
gsell's avatar
gsell committed
544 545
     */
    int  fieldflag;
546
    if(type == std::string("CARBONCYCL")) {
gsell's avatar
gsell committed
547
        fieldflag = 2;
548
    } else if(type == std::string("CYCIAE")) {
gsell's avatar
gsell committed
549
        fieldflag = 3;
550
    } else if(type == std::string("AVFEQ")) {
gsell's avatar
gsell committed
551
        fieldflag = 4;
552
    } else if(type == std::string("FFAG")) {
gsell's avatar
gsell committed
553
        fieldflag = 5;
554
    } else if(type == std::string("BANDRF")) {
gsell's avatar
gsell committed
555
        fieldflag = 6;
556 557
    } else if(type == std::string("SYNCHROCYCLOTRON")) {
	fieldflag = 7;
558
    } else //(type == "RING")
gsell's avatar
gsell committed
559 560
        fieldflag = 1;

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

    double BcParameter[8];
565

566
    for(int i = 0; i < 8; i++)
567
        BcParameter[i] = 0.0;
568

569 570
    BcParameter[0] = 0.001 * elptr->getRmin();
    BcParameter[1] = 0.001 * elptr->getRmax();
gsell's avatar
gsell committed
571

572
    // Store inner radius and outer radius of cyclotron field map in the list
573
    buildupFieldList(BcParameter, ElementBase::CYCLOTRON, elptr);
gsell's avatar
gsell committed
574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590
}

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

/**
 *
 *
 * @param coll
 */
void ParallelCyclotronTracker::visitCollimator(const Collimator &coll) {

591
    *gmsg << "* --------- Collimator -----------------------------" << endl;
gsell's avatar
gsell committed
592

593 594
    Collimator* elptr = dynamic_cast<Collimator *>(coll.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
595

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

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

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

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

608
    double zstart = elptr->getZStart();
adelmann's avatar
adelmann committed
609
    *gmsg << "* Zstart= " << zstart << " [mm]" << endl;
610 611

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

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

frey_m's avatar
frey_m committed
617
    elptr->initialise(itsBunch_m);
gsell's avatar
gsell committed
618 619 620 621

    double BcParameter[8];
    for(int i = 0; i < 8; i++)
        BcParameter[i] = 0.0;
622

623 624 625 626 627 628
    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 ;

629
    buildupFieldList(BcParameter, ElementBase::COLLIMATOR, elptr);
gsell's avatar
gsell committed
630 631 632 633 634 635 636 637 638 639 640 641
}

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

/**
 *
 *
 * @param lamb
 */
void ParallelCyclotronTracker::visitLambertson(const Lambertson &lamb) {
    *gmsg << "In Lambertson; L= " << lamb.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Lambertson *>(lamb.clone()));
}

684 685 686
void ParallelCyclotronTracker::visitOffset(const Offset & off) {
    if (opalRing_m == NULL)
        throw OpalException(
687 688
                            "ParallelCylcotronTracker::visitOffset",
                            "Attempt to place an offset when Ring not defined");
689 690
    Offset* offNonConst = const_cast<Offset*>(&off);
    offNonConst->updateGeometry(opalRing_m->getNextPosition(),
691
                                opalRing_m->getNextNormal());
692 693 694
    opalRing_m->appendElement(off);
}

gsell's avatar
gsell committed
695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716
/**
 *
 *
 * @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());
}

717

gsell's avatar
gsell committed
718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733
/**
 *
 *
 * @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()));
}

/**
 *
 *
 * @param prob
 */
void ParallelCyclotronTracker::visitProbe(const Probe &prob) {
734
    *gmsg << "* -----------  Probe -------------------------------" << endl;
735 736
    Probe *elptr = dynamic_cast<Probe *>(prob.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
737

738
    double xstart = elptr->getXstart();
739
    *gmsg << "XStart= " << xstart << " [mm]" << endl;
gsell's avatar
gsell committed
740

741
    double xend = elptr->getXend();
742
    *gmsg << "XEnd= " << xend << " [mm]" << endl;
gsell's avatar
gsell committed
743

744
    double ystart = elptr->getYstart();
745
    *gmsg << "YStart= " << ystart << " [mm]" << endl;
gsell's avatar
gsell committed
746

747
    double yend = elptr->getYend();
748
    *gmsg << "YEnd= " << yend << " [mm]" << endl;
gsell's avatar
gsell committed
749

750
    double width = elptr->getWidth();
751
    *gmsg << "Width= " << width << " [mm]" << endl;
gsell's avatar
gsell committed
752 753 754


    // initialise, do nothing
frey_m's avatar
frey_m committed
755
    elptr->initialise(itsBunch_m);
gsell's avatar
gsell committed
756 757 758 759

    double BcParameter[8];
    for(int i = 0; i < 8; i++)
        BcParameter[i] = 0.0;
760

761 762 763 764 765
    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 ;
gsell's avatar
gsell committed
766 767

    // store probe parameters in the list
768
    buildupFieldList(BcParameter, ElementBase::PROBE, elptr);
gsell's avatar
gsell committed
769 770 771 772 773 774 775 776 777 778 779 780
}

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

781 782 783 784 785 786
void ParallelCyclotronTracker::visitSBend3D(const SBend3D &bend) {
    *gmsg << "Adding SBend3D" << endl;
    if (opalRing_m != NULL)
        opalRing_m->appendElement(bend);
    else
        throw OpalException("ParallelCyclotronTracker::visitSBend3D",
787
                            "Need to define a RINGDEFINITION to use SBend3D element");
788 789
}

790 791
void ParallelCyclotronTracker::visitScalingFFAGMagnet(const ScalingFFAGMagnet &bend) {
    *gmsg << "Adding ScalingFFAGMagnet" << endl;
ext-rogers_c's avatar
ext-rogers_c committed
792
    if (opalRing_m != NULL) {
793
        opalRing_m->appendElement(bend);
ext-rogers_c's avatar
ext-rogers_c committed
794
    } else {
795 796
        throw OpalException("ParallelCyclotronTracker::visitScalingFFAGMagnet",
                            "Need to define a RINGDEFINITION to use ScalingFFAGMagnet element");
ext-rogers_c's avatar
ext-rogers_c committed
797
    }
798 799
}

800 801 802 803 804 805
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",
806
                            "Need to define a RINGDEFINITION to use VariableRFCavity element");
807 808
}

gsell's avatar
gsell committed
809 810 811 812 813 814 815 816
/**
 *
 *
 * @param as
 */
void ParallelCyclotronTracker::visitRFCavity(const RFCavity &as) {

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

818 819
    RFCavity *elptr = dynamic_cast<RFCavity *>(as.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
820 821 822 823 824 825 826

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

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

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

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

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

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

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

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

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

851 852 853
    /*
      Setup time dependence and in case of no
      timedependence use a polynom with  a_0 = 1 and a_k = 0, k = 1,2,3.
854
    */
855

856 857 858 859
    std::shared_ptr<AbstractTimeDependence> freq_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> ampl_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> phase_atd = nullptr;

860
    dvector_t  unityVec;
861 862 863 864
    unityVec.push_back(1.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
865

866
    if (elptr->getFrequencyModelName() != "") {
867 868
        freq_atd = AbstractTimeDependence::getTimeDependence(elptr->getFrequencyModelName());
        *gmsg << "* Variable frequency RF Model name " << elptr->getFrequencyModelName() << endl;
869
    }
870 871
    else
        freq_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
872 873

    if (elptr->getAmplitudeModelName() != "") {
874 875
        ampl_atd = AbstractTimeDependence::getTimeDependence(elptr->getAmplitudeModelName());
        *gmsg << "* Variable amplitude RF Model name " << elptr->getAmplitudeModelName() << endl;
876
    }
877 878
    else
        ampl_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
879 880

    if (elptr->getPhaseModelName() != "") {
881 882
        phase_atd = AbstractTimeDependence::getTimeDependence(elptr->getPhaseModelName());
        *gmsg << "* Variable phase RF Model name " << elptr->getPhaseModelName() << endl;
883
    }
884 885
    else
        phase_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
886

gsell's avatar
gsell committed
887
    // read cavity voltage profile data from file.
frey_m's avatar
frey_m committed
888
    elptr->initialise(itsBunch_m, freq_atd, ampl_atd, phase_atd);
gsell's avatar
gsell committed
889 890 891 892

    double BcParameter[8];
    for(int i = 0; i < 8; i++)
        BcParameter[i] = 0.0;
893

894 895
    BcParameter[0] = 0.001 * rmin;
    BcParameter[1] = 0.001 * rmax;
896
    BcParameter[2] = 0.001 * pdis;
gsell's avatar
gsell committed
897 898
    BcParameter[3] = angle;

899
    buildupFieldList(BcParameter, ElementBase::RFCAVITY, elptr);
gsell's avatar
gsell committed
900 901 902 903 904 905 906 907
}

/**
 *
 *
 * @param rfq
 */
void ParallelCyclotronTracker::visitRFQuadrupole(const RFQuadrupole &rfq) {
908
    *gmsg << "In RFQuadrupole; L = " << rfq.getElementLength() << " however the element is missing " << endl;
gsell's avatar
gsell committed
909 910 911 912 913 914 915 916 917
    myElements.push_back(dynamic_cast<RFQuadrupole *>(rfq.clone()));
}

/**
 *
 *
 * @param bend
 */
void ParallelCyclotronTracker::visitSBend(const SBend &bend) {