ParallelCyclotronTracker.cpp 228 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 98
const double c_mmtns = Physics::c * 1.0e-6; // m/s --> mm/ns
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
99

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

#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,
119 120
                                                   const PartData &reference,
                                                   bool revBeam, bool revTrack):
gsell's avatar
gsell committed
121
    Tracker(beamline, reference, revBeam, revTrack),
122
    eta_m(0.01),
gsell's avatar
gsell committed
123 124
    myNode_m(Ippl::myNode()),
    initialLocalNum_m(0),
125
    initialTotalNum_m(0),
126 127
    opalRing_m(NULL),
    lastDumpedStep_m(0) {
gsell's avatar
gsell committed
128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147
    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,
148
                                                   int maxSTEPS, int timeIntegrator):
gsell's avatar
gsell committed
149 150 151
    Tracker(beamline, reference, revBeam, revTrack),
    maxSteps_m(maxSTEPS),
    timeIntegrator_m(timeIntegrator),
152
    eta_m(0.01),
gsell's avatar
gsell committed
153 154
    myNode_m(Ippl::myNode()),
    initialLocalNum_m(bunch.getLocalNum()),
155
    initialTotalNum_m(bunch.getTotalNum()),
156
    opalRing_m(NULL),
157
    lastDumpedStep_m(0) {
gsell's avatar
gsell committed
158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175
    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() {
176
    for(std::list<Component *>::iterator compindex = myElements.begin(); compindex != myElements.end(); compindex++) {
gsell's avatar
gsell committed
177 178 179 180 181 182
        delete(*compindex);
    }
    for(beamline_list::iterator fdindex = FieldDimensions.begin(); fdindex != FieldDimensions.end(); fdindex++) {
        delete(*fdindex);
    }
    delete itsBeamline;
183 184 185
    if (opalRing_m != NULL) {
        // delete opalRing_m;
    }
gsell's avatar
gsell committed
186 187
}

188 189 190 191 192 193
/**
 * AAA
 *
 * @param none
 */
void ParallelCyclotronTracker::initializeBoundaryGeometry() {
194 195 196 197 198 199 200 201 202 203 204 205
    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;
    }
206
}
207

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

216
    Inform msg("bgf_main_collision_test ");
217

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

222
    Vector_t intecoords = 0.0;
223

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

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


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

245
    std::string  SfileName2 = SfileName + std::string("-Angle0.dat");
gsell's avatar
gsell committed
246 247

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    if(spiral_flag) {

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

    }

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

386 387 388 389 390
        // Get reference values from cyclotron element
        referenceR     = elptr->getRinit();
        referenceTheta = elptr->getPHIinit();
        referenceZ     = elptr->getZinit();
        referencePr    = elptr->getPRinit();
391 392
        referencePz    = elptr->getPZinit();

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

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

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

        if(insqrt < 0) {

            if(insqrt > -1.0e-10) {

	        referencePt = 0.0;

            } else {

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

        } else {

            referencePt = sqrt(insqrt);
        }

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

425
        // Restart a run:
426 427
    } else {

428
        // If the user wants to save the restarted run in local frame,
429 430
        // make sure the previous h5 file was local too
        if (Options::psDumpLocalFrame) {
431

432
	    if (!previousH5Local) {
433

434 435
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                    "You are trying a local restart from a global h5 file!");
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
439
        } else {
440

441
	    if (previousH5Local) {
442

443 444
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                    "You are trying a global restart from a local h5 file!");
445
            }
446
        }
447

448
        // Adjust some of the reference variables from the h5 file
449 450
        referencePhi *= Physics::deg2rad;
        referencePsi *= Physics::deg2rad;
451
        referencePtot = bega;
452 453 454

        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {

455
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
456 457
                                "PHIINIT is out of [-180, 180)!");
        }
458 459
    }

460 461
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
462

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

gsell's avatar
gsell committed
478
    double sym = elptr->getSymmetry();
479
    *gmsg << "* " << sym << "-fold field symmerty " << endl;
gsell's avatar
gsell committed
480

481 482 483
    // 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
484

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

488
    std::string type = elptr->getCyclotronType();
489
    *gmsg << "* Type of cyclotron = " << type << " " << endl;
490

491 492
    double rmin = elptr->getMinR();
    double rmax = elptr->getMaxR();
493
    *gmsg << "* Radial aperture = " << rmin << " ... " << rmax<<" [mm] "<< endl;
494 495 496

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

499
    /**
500 501 502 503 504 505 506 507
       bool Sflag = elptr->getSuperpose();
       std::string flagsuperposed;
       if (Sflag)
       flagsuperposed="yes";
       else
       flagsuperposed="no";
       *gmsg << "* Electric field maps are superposed? " << flagsuperposed << " " << endl;
       */
508

gsell's avatar
gsell committed
509
    double h = elptr->getCyclHarm();
510
    *gmsg << "* Number of trimcoils = " << elptr->getNumberOfTrimcoils() << endl;
511
    *gmsg << "* Harmonic number h = " << h << " " << endl;
512

513
    /**
514 515 516
       if (elptr->getSuperpose())
       *gmsg << "* Fields are superposed " << endl;
       */
adelmann's avatar
adelmann committed
517

gsell's avatar
gsell committed
518 519 520 521 522 523 524 525 526 527 528
    /**
     * 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
     */
    int  fieldflag;
529
    if(type == std::string("CARBONCYCL")) {
gsell's avatar
gsell committed
530
        fieldflag = 2;
531
    } else if(type == std::string("CYCIAE")) {
gsell's avatar
gsell committed
532
        fieldflag = 3;
533
    } else if(type == std::string("AVFEQ")) {
gsell's avatar
gsell committed
534
        fieldflag = 4;
535
    } else if(type == std::string("FFAG")) {
gsell's avatar
gsell committed
536
        fieldflag = 5;
537
    } else if(type == std::string("BANDRF")) {
gsell's avatar
gsell committed
538
        fieldflag = 6;
539 540
    } else if(type == std::string("SYNCHROCYCLOTRON")) {
	fieldflag = 7;
541
    } else //(type == "RING")
gsell's avatar
gsell committed
542 543 544 545
        fieldflag = 1;

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

    double BcParameter[8];
550

551
    for(int i = 0; i < 8; i++)
552
        BcParameter[i] = 0.0;
553

gsell's avatar
gsell committed
554 555 556
    BcParameter[0] = elptr->getRmin();
    BcParameter[1] = elptr->getRmax();

557
    // Store inner radius and outer radius of cyclotron field map in the list
558
    buildupFieldList(BcParameter, ElementBase::CYCLOTRON, elptr);
gsell's avatar
gsell committed
559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575
}

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

576
    *gmsg << "* --------- Collimator -----------------------------" << endl;
gsell's avatar
gsell committed
577

578 579
    Collimator* elptr = dynamic_cast<Collimator *>(coll.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
580

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

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

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

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

593
    double zstart = elptr->getZStart();
adelmann's avatar
adelmann committed
594
    *gmsg << "* Zstart= " << zstart << " [mm]" << endl;
595 596

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

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

602
    elptr->initialise(itsBunch);
gsell's avatar
gsell committed
603 604 605 606

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

608 609 610 611
    BcParameter[0] = xstart ;
    BcParameter[1] = xend;
    BcParameter[2] = ystart ;
    BcParameter[3] = yend;
gsell's avatar
gsell committed
612
    BcParameter[4] = width ;
613
    buildupFieldList(BcParameter, ElementBase::COLLIMATOR, elptr);
gsell's avatar
gsell committed
614 615 616 617 618 619 620 621 622 623 624 625
}

/**
 *
 *
 * @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
626 627 628 629 630 631 632 633 634 635 636 637
/**
 *
 *
 * @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
638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667
/**
 *
 *
 * @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()));
}

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

gsell's avatar
gsell committed
679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700
/**
 *
 *
 * @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());
}

701

gsell's avatar
gsell committed
702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717
/**
 *
 *
 * @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) {
718
    *gmsg << "* -----------  Probe -------------------------------" << endl;
719 720
    Probe *elptr = dynamic_cast<Probe *>(prob.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
721

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

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

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

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

734
    double width = elptr->getWidth();
735
    *gmsg << "Width= " << width << " [mm]" << endl;
gsell's avatar
gsell committed
736 737 738


    // initialise, do nothing
739
    elptr->initialise(itsBunch);
gsell's avatar
gsell committed
740 741 742 743

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

gsell's avatar
gsell committed
745 746 747 748 749 750 751
    BcParameter[0] = xstart ;
    BcParameter[1] = xend;
    BcParameter[2] = ystart ;
    BcParameter[3] = yend;
    BcParameter[4] = width ;

    // store probe parameters in the list
752
    buildupFieldList(BcParameter, ElementBase::PROBE, elptr);
gsell's avatar
gsell committed
753 754 755 756 757 758 759 760 761 762 763 764
}

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

765 766 767 768 769 770
void ParallelCyclotronTracker::visitSBend3D(const SBend3D &bend) {
    *gmsg << "Adding SBend3D" << endl;
    if (opalRing_m != NULL)
        opalRing_m->appendElement(bend);
    else
        throw OpalException("ParallelCyclotronTracker::visitSBend3D",
771
                            "Need to define a RINGDEFINITION to use SBend3D element");
772 773
}

774 775 776 777 778 779
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",
780
                            "Need to define a RINGDEFINITION to use VariableRFCavity element");
781 782
}

gsell's avatar
gsell committed
783 784 785 786 787 788 789 790
/**
 *
 *
 * @param as
 */
void ParallelCyclotronTracker::visitRFCavity(const RFCavity &as) {

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

792 793
    RFCavity *elptr = dynamic_cast<RFCavity *>(as.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
794 795 796 797 798 799 800

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

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

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

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

810
    std::string fmfn = elptr->getFieldMapFN();
gsell's avatar
gsell committed
811
    *gmsg << "* RF Field map file name= " << fmfn << endl;
812
    double angle = elptr->getAzimuth();
gsell's avatar
gsell committed
813 814
    *gmsg << "* Cavity azimuth position= " << angle << " [deg] " << endl;

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

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

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

824 825 826 827

    /*
      Setup time dependence and in case of no
      timedependence use a polynom with  a_0 = 1 and a_k = 0, k = 1,2,3.
828
    */
829

830 831 832 833
    std::shared_ptr<AbstractTimeDependence> freq_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> ampl_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> phase_atd = nullptr;

834 835 836 837 838
    std::vector<double>  unityVec;
    unityVec.push_back(1.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
839

840
    if (elptr->getFrequencyModelName() != "") {
841 842
        freq_atd = AbstractTimeDependence::getTimeDependence(elptr->getFrequencyModelName());
        *gmsg << "* Variable frequency RF Model name " << elptr->getFrequencyModelName() << endl;
843
    }
844 845
    else
        freq_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
846 847

    if (elptr->getAmplitudeModelName() != "") {
848 849
        ampl_atd = AbstractTimeDependence::getTimeDependence(elptr->getAmplitudeModelName());
        *gmsg << "* Variable amplitude RF Model name " << elptr->getAmplitudeModelName() << endl;
850
    }
851 852
    else
        ampl_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
853 854

    if (elptr->getPhaseModelName() != "") {
855 856
        phase_atd = AbstractTimeDependence::getTimeDependence(elptr->getPhaseModelName());
        *gmsg << "* Variable phase RF Model name " << elptr->getPhaseModelName() << endl;
857
    }
858 859
    else
        phase_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
860

gsell's avatar
gsell committed
861
    // read cavity voltage profile data from file.
862
    elptr->initialise(itsBunch, freq_atd, ampl_atd, phase_atd);
gsell's avatar
gsell committed
863 864 865 866

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

gsell's avatar
gsell committed
868 869 870 871 872
    BcParameter[0] = rmin;
    BcParameter[1] = rmax;
    BcParameter[2] = pdis;
    BcParameter[3] = angle;

873
    buildupFieldList(BcParameter, ElementBase::RFCAVITY, elptr);
gsell's avatar
gsell committed
874 875 876 877 878 879 880 881
}

/**
 *
 *
 * @param rfq
 */
void ParallelCyclotronTracker::visitRFQuadrupole(const RFQuadrupole &rfq) {
882
    *gmsg << "In RFQuadrupole; L = " << rfq.getElementLength() << " however the element is missing " << endl;
gsell's avatar
gsell committed
883 884 885 886 887 888 889 890 891
    myElements.push_back(dynamic_cast<RFQuadrupole *>(rfq.clone()));
}

/**
 *
 *
 * @param bend
 */
void ParallelCyclotronTracker::visitSBend(const SBend &bend) {
892
    *gmsg << "In SBend; L = " << bend.getElementLength() << " however the element is missing " << endl;
gsell's avatar
gsell committed
893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911
    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) {
912

913
    *gmsg << endl << "* -----------------------------  Septum ------------------------------- *" << endl;
gsell's avatar
gsell committed
914

915 916
    Septum *elptr = dynamic_cast<Septum *>(sept.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
917

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

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

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

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

930
    double width = elptr->getWidth();
931
    *gmsg << "Width = " << width << " [mm]" << endl;
gsell's avatar
gsell committed
932 933 934


    // initialise, do nothing
935
    elptr->initialise(itsBunch);
gsell's avatar
gsell committed
936 937 938 939

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

gsell's avatar
gsell committed
941 942 943 944 945 946 947
    BcParameter[0] = xstart ;
    BcParameter[1] = xend;
    BcParameter[2] = ystart ;
    BcParameter[3] = yend;
    BcParameter[4] = width ;

    // store septum parameters in the list