ParallelCyclotronTracker.cpp 154 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"
gsell's avatar
gsell committed
48 49 50 51 52
#include "AbsBeamline/Separator.h"
#include "AbsBeamline/Septum.h"
#include "AbsBeamline/Solenoid.h"
#include "AbsBeamline/CyclotronValley.h"
#include "AbsBeamline/Stripper.h"
53
#include "AbsBeamline/VariableRFCavity.h"
54

55 56
#include "AbstractObjects/Element.h"

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

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

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

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

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

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

class Beamline;
class PartData;
90

gsell's avatar
gsell committed
91 92 93 94 95 96
using Physics::m_p; // GeV
using Physics::PMASS;
using Physics::PCHARGE;
using Physics::pi;
using Physics::q_e;

97
const double c_mmtns = Physics::c * 1.0e-6; // m/s --> mm/ns
98
const double c_mtns = Physics::c * 1.0e-3;  // m/s --> m/ns
99
const double mass_coeff = 1.0e18 * q_e / Physics::c / Physics::c; // from GeV/c^2 to basic unit: GV*C*s^2/m^2
100

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

#define PSdim 6

extern Inform *gmsg;

// typedef FVector<double, PSdim> Vector;

/**
 * Constructor ParallelCyclotronTracker
 *
 * @param beamline
 * @param reference
 * @param revBeam
 * @param revTrack
 */
ParallelCyclotronTracker::ParallelCyclotronTracker(const Beamline &beamline,
120 121
                                                   const PartData &reference,
                                                   bool revBeam, bool revTrack):
gsell's avatar
gsell committed
122
    Tracker(beamline, reference, revBeam, revTrack),
123
    lastDumpedStep_m(0),
124
    eta_m(0.01),
gsell's avatar
gsell committed
125 126
    myNode_m(Ippl::myNode()),
    initialLocalNum_m(0),
127
    initialTotalNum_m(0),
128
    opalRing_m(NULL) {
gsell's avatar
gsell committed
129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148
    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,
                                                   PartBunch &bunch,
                                                   DataSink &ds,
                                                   const PartData &reference,
                                                   bool revBeam, bool revTrack,
149
                                                   int maxSTEPS, int timeIntegrator):
gsell's avatar
gsell committed
150 151
    Tracker(beamline, reference, revBeam, revTrack),
    maxSteps_m(maxSTEPS),
152
    lastDumpedStep_m(0),
gsell's avatar
gsell committed
153
    timeIntegrator_m(timeIntegrator),
154
    eta_m(0.01),
gsell's avatar
gsell committed
155 156
    myNode_m(Ippl::myNode()),
    initialLocalNum_m(bunch.getLocalNum()),
157
    initialTotalNum_m(bunch.getTotalNum()),
158
    opalRing_m(NULL) {
gsell's avatar
gsell committed
159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176
    itsBeamline = dynamic_cast<Beamline *>(beamline.clone());
    itsBunch = &bunch;
    itsDataSink = &ds;
    //  scaleFactor_m = itsBunch->getdT() * c;
    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");
}

/**
 * Destructor ParallelCyclotronTracker
 *
 */
ParallelCyclotronTracker::~ParallelCyclotronTracker() {
177
    for(std::list<Component *>::iterator compindex = myElements.begin(); compindex != myElements.end(); compindex++) {
gsell's avatar
gsell committed
178 179 180 181 182 183
        delete(*compindex);
    }
    for(beamline_list::iterator fdindex = FieldDimensions.begin(); fdindex != FieldDimensions.end(); fdindex++) {
        delete(*fdindex);
    }
    delete itsBeamline;
184 185 186
    if (opalRing_m != NULL) {
        // delete opalRing_m;
    }
gsell's avatar
gsell committed
187 188
}

189 190 191 192 193 194
/**
 * AAA
 *
 * @param none
 */
void ParallelCyclotronTracker::initializeBoundaryGeometry() {
195 196 197 198 199 200 201 202 203 204 205 206
    for(std::list<Component *>::iterator compindex = myElements.begin(); compindex != myElements.end(); compindex++) {
        bgf_m = dynamic_cast<ElementBase *>(*compindex)->getBoundaryGeometry();
        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;
    }
207
}
208

209 210 211 212 213 214
/**
 *
 *
 * @param fn Base file name
 */
void ParallelCyclotronTracker::bgf_main_collision_test() {
215
    if(!bgf_m) return;
216

217
    Inform msg("bgf_main_collision_test ");
218

219 220 221
    /**
     *Here we check if a particles is outside the domain, flag it for deletion
     */
222

223
    Vector_t intecoords = 0.0;
224

225 226 227
    // This has to match the dT in the rk4 pusher! -DW
    //double dtime = 0.5 * itsBunch->getdT();  // Old
    double dtime = itsBunch->getdT() * getHarmonicNumber();  // New
228

229 230
    int triId = 0;
    for(size_t i = 0; i < itsBunch->getLocalNum(); i++) {
231 232
        int res = bgf_m->partInside(itsBunch->R[i], itsBunch->P[i], dtime, itsBunch->PType[i], itsBunch->Q[i], intecoords, triId);
        //int res = bgf_m->partInside(itsBunch->R[i]*1.0e-3, itsBunch->P[i], dtime, itsBunch->PType[i], itsBunch->Q[i], intecoords, triId);
233 234 235
        if(res >= 0) {
            itsBunch->Bin[i] = -1;
        }
236
    }
237 238 239
}


gsell's avatar
gsell committed
240 241 242 243 244
/**
 *
 *
 * @param fn Base file name
 */
245
void ParallelCyclotronTracker::openFiles(std::string SfileName) {
gsell's avatar
gsell committed
246

247
    std::string  SfileName2 = SfileName + std::string("-Angle0.dat");
gsell's avatar
gsell committed
248
    outfTheta0_m.precision(8);
249
    outfTheta0_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
250
    outfTheta0_m.open(SfileName2.c_str());
251
    outfTheta0_m << "#  r [mm]      beta_r*gamma       theta [mm]      beta_theta*gamma        z [mm]          beta_z*gamma" << std::endl;
gsell's avatar
gsell committed
252

253
    SfileName2 = SfileName + std::string("-Angle1.dat");
gsell's avatar
gsell committed
254
    outfTheta1_m.precision(8);
255
    outfTheta1_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
256
    outfTheta1_m.open(SfileName2.c_str());
257
    outfTheta1_m << "#  r [mm]      beta_r*gamma       theta [mm]      beta_theta*gamma        z [mm]          beta_z*gamma"  << std::endl;
gsell's avatar
gsell committed
258

259
    SfileName2 = SfileName + std::string("-Angle2.dat");
gsell's avatar
gsell committed
260
    outfTheta2_m.precision(8);
261
    outfTheta2_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
262
    outfTheta2_m.open(SfileName2.c_str());
263
    outfTheta2_m << "#  r [mm]      beta_r*gamma       theta [mm]      beta_theta*gamma        z [mm]          beta_z*gamma"  << std::endl;
gsell's avatar
gsell committed
264 265 266

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

267
    SfileName2 = SfileName + std::string("-afterEachTurn.dat");
gsell's avatar
gsell committed
268 269

    outfThetaEachTurn_m.precision(8);
270
    outfThetaEachTurn_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
271 272

    outfThetaEachTurn_m.open(SfileName2.c_str());
273
    outfTheta2_m << "#  r [mm]      beta_r*gamma       theta [mm]      beta_theta*gamma        z [mm]          beta_z*gamma"  << std::endl;
gsell's avatar
gsell committed
274 275 276 277 278 279 280 281 282 283 284 285 286 287 288
}

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

289
/**
290 291 292
 *
 * @param ring
 */
kraus's avatar
kraus committed
293
void ParallelCyclotronTracker::visitRing(const Ring &ring) {
Daniel Winklehner's avatar
Daniel Winklehner committed
294

295
    *gmsg << "* ----------------------------- Adding Ring ------------------------------ *" << endl;
Daniel Winklehner's avatar
Daniel Winklehner committed
296 297

    if (opalRing_m != NULL)
298
        delete opalRing_m;
Daniel Winklehner's avatar
Daniel Winklehner committed
299

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

302
    myElements.push_back(opalRing_m);
Daniel Winklehner's avatar
Daniel Winklehner committed
303

304 305 306 307 308
    opalRing_m->initialise(itsBunch);

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

310
    if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
kraus's avatar
kraus committed
311
        throw OpalException("Error in ParallelCyclotronTracker::visitRing",
312 313
                            "PHIINIT is out of [-180, 180)!");
    }
Daniel Winklehner's avatar
Daniel Winklehner committed
314 315

    referenceZ = 0.0;
316
    referencePz = 0.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
317 318 319 320

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

321 322
    if(referencePtot < 0.0)
        referencePt *= -1.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
323

324 325
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
gsell's avatar
gsell committed
326

327
    double BcParameter[8];
Daniel Winklehner's avatar
Daniel Winklehner committed
328

329
    for(int i = 0; i < 8; i++)
Daniel Winklehner's avatar
Daniel Winklehner committed
330 331
        BcParameter[i] = 0.0;

332
    buildupFieldList(BcParameter, ElementBase::RING, opalRing_m);
333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348

    // 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
349 350 351 352 353 354 355 356

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

357
    *gmsg << "* -------------------------- Adding Cyclotron ---------------------------- *" << endl;
gsell's avatar
gsell committed
358

359 360
    Cyclotron *elptr = dynamic_cast<Cyclotron *>(cycl.clone());
    myElements.push_back(elptr);
361

362
    // Is this a Spiral Inflector Simulation? If yes, we'll give the user some
363 364 365 366 367 368
    // useful information
    spiral_flag = elptr->getSpiralFlag();

    if(spiral_flag) {

        *gmsg << endl << "* This is a Spiral Inflector Simulation! This means the following:" << endl;
369
        *gmsg         << "* 1.) It is up to the user to provide appropriate geometry, electric and magnetic fields!" << endl;
370 371
        *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;
372
        *gmsg         << "* 2.) For high currentst is strongly recommended to use the SAAMG fieldsolver," << endl;
373 374
        *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;
375 376
        *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;
377 378 379 380 381 382 383
        *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;

    }

384
    // Fresh run (no restart):
385
    if(!OpalData::getInstance()->inRestartRun()) {
386

387 388
        // Get reference values from cyclotron element 
        // For now, these are still stored in mm. should be the only ones. -DW
389 390 391 392
        referenceR     = elptr->getRinit();
        referenceTheta = elptr->getPHIinit();
        referenceZ     = elptr->getZinit();
        referencePr    = elptr->getPRinit();
393 394
        referencePz    = elptr->getPZinit();

395
        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
396 397
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                "PHIINIT is out of [-180, 180)!");
398 399 400 401 402
        }

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

        // Calculate reference azimuthal (tangential) momentum from total-, z- and radial momentum:
403
        float insqrt = referencePtot * referencePtot - \
404
            referencePr * referencePr - referencePz * referencePz;
405 406 407 408 409 410 411 412 413

        if(insqrt < 0) {

            if(insqrt > -1.0e-10) {

	        referencePt = 0.0;

            } else {

414
	        throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
415
                                    "Pt imaginary!");
416 417 418 419 420 421 422
            }

        } else {

            referencePt = sqrt(insqrt);
        }

423
        if(referencePtot < 0.0)
Daniel Winklehner's avatar
Daniel Winklehner committed
424
            referencePt *= -1.0;
425 426
        // End calculate referencePt

427
        // Restart a run:
428 429
    } else {

430
        // If the user wants to save the restarted run in local frame,
431
        // make sure the previous h5 file was local too
winklehner_d's avatar
-DW  
winklehner_d committed
432 433
      if (Options::psDumpLocalFrame != Options::GLOBAL) {
        if (!previousH5Local) {
434 435
                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
436
        }
437 438
            // 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
439 440
      } else {
        if (previousH5Local) {
441 442
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                    "You are trying a global restart from a local h5 file!");
443
            }
444
        }
445

446
        // Adjust some of the reference variables from the h5 file
447 448
        referencePhi *= Physics::deg2rad;
        referencePsi *= Physics::deg2rad;
449
        referencePtot = bega;
450
        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
451
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
452 453
                                "PHIINIT is out of [-180, 180)!");
        }
454 455
    }

456 457
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
458

459
    *gmsg << endl;
adelmann's avatar
adelmann committed
460
    *gmsg << "* Bunch global starting position:" << endl;
461 462
    *gmsg << "* RINIT = " << referenceR  << " [mm]" << endl;
    *gmsg << "* PHIINIT = " << referenceTheta << " [deg]" << endl;
463
    *gmsg << "* ZINIT = " << referenceZ << " [mm]" << endl;
464
    *gmsg << endl;
adelmann's avatar
adelmann committed
465
    *gmsg << "* Bunch global starting momenta:" << endl;
gsell's avatar
gsell committed
466 467
    *gmsg << "* Initial gamma = " << itsReference.getGamma() << endl;
    *gmsg << "* Initial beta = " << itsReference.getBeta() << endl;
468
    *gmsg << "* Reference total momentum (beta * gamma) = " << referencePtot * 1000.0 << " [MCU]" << endl;
469 470 471
    *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;
472
    *gmsg << endl;
adelmann's avatar
adelmann committed
473

gsell's avatar
gsell committed
474
    double sym = elptr->getSymmetry();
475
    *gmsg << "* " << sym << "-fold field symmerty " << endl;
gsell's avatar
gsell committed
476

477 478 479
    // 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
480

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

484
    std::string type = elptr->getCyclotronType();
485
    *gmsg << "* Type of cyclotron = " << type << " " << endl;
486

487 488
    double rmin = elptr->getMinR();
    double rmax = elptr->getMaxR();
489
    *gmsg << "* Radial aperture = " << rmin << " ... " << rmax<<" [mm] "<< endl;
490 491 492

    double zmin = elptr->getMinZ();
    double zmax = elptr->getMaxZ();
493
    *gmsg << "* Vertical aperture = " << zmin << " ... " << zmax<<" [mm]"<< endl;
gsell's avatar
gsell committed
494 495

    double h = elptr->getCyclHarm();
496
    *gmsg << "* Number of trimcoils = " << elptr->getNumberOfTrimcoils() << endl;
497
    *gmsg << "* Harmonic number h = " << h << " " << endl;
498

gsell's avatar
gsell committed
499 500 501 502 503 504 505 506 507
    /**
     * 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
508
     * fieldflag = 7, read in fields for Daniel's synchrocyclotron simulations
gsell's avatar
gsell committed
509 510
     */
    int  fieldflag;
511
    if(type == std::string("CARBONCYCL")) {
gsell's avatar
gsell committed
512
        fieldflag = 2;
513
    } else if(type == std::string("CYCIAE")) {
gsell's avatar
gsell committed
514
        fieldflag = 3;
515
    } else if(type == std::string("AVFEQ")) {
gsell's avatar
gsell committed
516
        fieldflag = 4;
517
    } else if(type == std::string("FFAG")) {
gsell's avatar
gsell committed
518
        fieldflag = 5;
519
    } else if(type == std::string("BANDRF")) {
gsell's avatar
gsell committed
520
        fieldflag = 6;
521 522
    } else if(type == std::string("SYNCHROCYCLOTRON")) {
	fieldflag = 7;
523
    } else //(type == "RING")
gsell's avatar
gsell committed
524 525 526 527
        fieldflag = 1;

    // read field map on the  middle plane of cyclotron.
    // currently scalefactor is set to 1.0
528
    // TEMP changed 1.0 to getBScale() to test if we can scale the midplane field -DW
529
    elptr->initialise(itsBunch, fieldflag, elptr->getBScale());
gsell's avatar
gsell committed
530 531

    double BcParameter[8];
532

533
    for(int i = 0; i < 8; i++)
534
        BcParameter[i] = 0.0;
535

536 537
    BcParameter[0] = 0.001 * elptr->getRmin();
    BcParameter[1] = 0.001 * elptr->getRmax();
gsell's avatar
gsell committed
538

539
    // Store inner radius and outer radius of cyclotron field map in the list
540
    buildupFieldList(BcParameter, ElementBase::CYCLOTRON, elptr);
gsell's avatar
gsell committed
541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557
}

/**
 * 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) {

558
    *gmsg << "* --------- Collimator -----------------------------" << endl;
gsell's avatar
gsell committed
559

560 561
    Collimator* elptr = dynamic_cast<Collimator *>(coll.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
562

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

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

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

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

575
    double zstart = elptr->getZStart();
adelmann's avatar
adelmann committed
576
    *gmsg << "* Zstart= " << zstart << " [mm]" << endl;
577 578

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

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

584
    elptr->initialise(itsBunch);
gsell's avatar
gsell committed
585 586 587 588

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

590 591 592 593 594 595
    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 ;

596
    buildupFieldList(BcParameter, ElementBase::COLLIMATOR, elptr);
gsell's avatar
gsell committed
597 598 599 600 601 602 603 604 605 606 607 608
}

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

651 652 653
void ParallelCyclotronTracker::visitOffset(const Offset & off) {
    if (opalRing_m == NULL)
        throw OpalException(
654 655
                            "ParallelCylcotronTracker::visitOffset",
                            "Attempt to place an offset when Ring not defined");
656 657
    Offset* offNonConst = const_cast<Offset*>(&off);
    offNonConst->updateGeometry(opalRing_m->getNextPosition(),
658
                                opalRing_m->getNextNormal());
659 660 661
    opalRing_m->appendElement(off);
}

gsell's avatar
gsell committed
662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683
/**
 *
 *
 * @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());
}

684

gsell's avatar
gsell committed
685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700
/**
 *
 *
 * @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) {
701
    *gmsg << "* -----------  Probe -------------------------------" << endl;
702 703
    Probe *elptr = dynamic_cast<Probe *>(prob.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
704

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

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

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

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

717
    double width = elptr->getWidth();
718
    *gmsg << "Width= " << width << " [mm]" << endl;
gsell's avatar
gsell committed
719 720 721


    // initialise, do nothing
722
    elptr->initialise(itsBunch);
gsell's avatar
gsell committed
723 724 725 726

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

728 729 730 731 732
    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
733 734

    // store probe parameters in the list
735
    buildupFieldList(BcParameter, ElementBase::PROBE, elptr);
gsell's avatar
gsell committed
736 737 738 739 740 741 742 743 744 745 746 747
}

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

748 749 750 751 752 753
void ParallelCyclotronTracker::visitSBend3D(const SBend3D &bend) {
    *gmsg << "Adding SBend3D" << endl;
    if (opalRing_m != NULL)
        opalRing_m->appendElement(bend);
    else
        throw OpalException("ParallelCyclotronTracker::visitSBend3D",
754
                            "Need to define a RINGDEFINITION to use SBend3D element");
755 756
}

757 758 759 760 761 762
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",
763
                            "Need to define a RINGDEFINITION to use VariableRFCavity element");
764 765
}

gsell's avatar
gsell committed
766 767 768 769 770 771 772 773
/**
 *
 *
 * @param as
 */
void ParallelCyclotronTracker::visitRFCavity(const RFCavity &as) {

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

775 776
    RFCavity *elptr = dynamic_cast<RFCavity *>(as.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
777 778 779 780 781 782 783

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

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

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

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

793
    std::string fmfn = elptr->getFieldMapFN();
gsell's avatar
gsell committed
794
    *gmsg << "* RF Field map file name= " << fmfn << endl;
795
    double angle = elptr->getAzimuth();
gsell's avatar
gsell committed
796 797
    *gmsg << "* Cavity azimuth position= " << angle << " [deg] " << endl;

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

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

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

807 808 809
    /*
      Setup time dependence and in case of no
      timedependence use a polynom with  a_0 = 1 and a_k = 0, k = 1,2,3.
810
    */
811

812 813 814 815
    std::shared_ptr<AbstractTimeDependence> freq_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> ampl_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> phase_atd = nullptr;

816 817 818 819 820
    std::vector<double>  unityVec;
    unityVec.push_back(1.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
821

822
    if (elptr->getFrequencyModelName() != "") {
823 824
        freq_atd = AbstractTimeDependence::getTimeDependence(elptr->getFrequencyModelName());
        *gmsg << "* Variable frequency RF Model name " << elptr->getFrequencyModelName() << endl;
825
    }
826 827
    else
        freq_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
828 829

    if (elptr->getAmplitudeModelName() != "") {
830 831
        ampl_atd = AbstractTimeDependence::getTimeDependence(elptr->getAmplitudeModelName());
        *gmsg << "* Variable amplitude RF Model name " << elptr->getAmplitudeModelName() << endl;
832
    }
833 834
    else
        ampl_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
835 836

    if (elptr->getPhaseModelName() != "") {
837 838
        phase_atd = AbstractTimeDependence::getTimeDependence(elptr->getPhaseModelName());
        *gmsg << "* Variable phase RF Model name " << elptr->getPhaseModelName() << endl;
839
    }
840 841
    else
        phase_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
842

gsell's avatar
gsell committed
843
    // read cavity voltage profile data from file.
844
    elptr->initialise(itsBunch, freq_atd, ampl_atd, phase_atd);
gsell's avatar
gsell committed
845 846 847 848

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

850 851
    BcParameter[0] = 0.001 * rmin;
    BcParameter[1] = 0.001 * rmax;
gsell's avatar
gsell committed
852 853 854
    BcParameter[2] = pdis;
    BcParameter[3] = angle;

855
    buildupFieldList(BcParameter, ElementBase::RFCAVITY, elptr);
gsell's avatar
gsell committed
856 857 858 859 860 861 862 863
}

/**
 *
 *
 * @param rfq
 */
void ParallelCyclotronTracker::visitRFQuadrupole(const RFQuadrupole &rfq) {
864
    *gmsg << "In RFQuadrupole; L = " << rfq.getElementLength() << " however the element is missing " << endl;
gsell's avatar
gsell committed
865 866 867 868 869 870 871 872 873
    myElements.push_back(dynamic_cast<RFQuadrupole *>(rfq.clone()));
}

/**
 *
 *
 * @param bend
 */
void ParallelCyclotronTracker::visitSBend(const SBend &bend) {
874
    *gmsg << "In SBend; L = " << bend.getElementLength() << " however the element is missing " << endl;
gsell's avatar
gsell committed
875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893
    myElements.push_back(dynamic_cast<SBend *>(bend.clone()));
}

/**
 *
 *
 * @param sep
 */
void ParallelCyclotronTracker::visitSeparator(const Separator &sep) {
    *gmsg << "In Seapator L= " << sep.getElementLength() << " however the element is missing " << endl;
    myElements.push_back(dynamic_cast<Separator *>(sep.clone()));
}

/**
 *
 *
 * @param sept
 */
void ParallelCyclotronTracker::visitSeptum(const Septum &sept) {
894

895
    *gmsg << endl << "* -----------------------------  Septum ------------------------------- *" << endl;
gsell's avatar
gsell committed
896

897 898
    Septum *elptr = dynamic_cast<Septum *>(sept.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
899

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

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

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

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

912
    double width = elptr->getWidth();
913
    *gmsg << "Width = " << width << " [mm]" << endl;
gsell's avatar
gsell committed
914 915 916


    // initialise, do nothing
917
    elptr->initialise(itsBunch);
gsell's avatar
gsell committed
918 919 920 921

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

923 924 925 926 927
    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
928 929

    // store septum parameters in the list
930
    buildupFieldList(BcParameter, ElementBase::SEPTUM, elptr);
gsell's avatar
gsell committed
931 932 933 934 935 936 937 938 939 940 941