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 432
        // make sure the previous h5 file was local too
        if (Options::psDumpLocalFrame) {
433

434
	    if (!previousH5Local) {
435

436 437
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                    "You are trying a local restart from a global h5 file!");
438
	    }
439 440
            // Else, if the user wants to save the restarted run in global frame,
            // make sure the previous h5 file was global too
441
        } else {
442

443
	    if (previousH5Local) {
444

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

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

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

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

462 463
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
464

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

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

483 484 485
    // 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
486

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

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

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

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

    double h = elptr->getCyclHarm();
502
    *gmsg << "* Number of trimcoils = " << elptr->getNumberOfTrimcoils() << endl;
503
    *gmsg << "* Harmonic number h = " << h << " " << endl;
504

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

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

    double BcParameter[8];
538

539
    for(int i = 0; i < 8; i++)
540
        BcParameter[i] = 0.0;
541

542 543
    BcParameter[0] = 0.001 * elptr->getRmin();
    BcParameter[1] = 0.001 * elptr->getRmax();
gsell's avatar
gsell committed
544

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

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

564
    *gmsg << "* --------- Collimator -----------------------------" << endl;
gsell's avatar
gsell committed
565

566 567
    Collimator* elptr = dynamic_cast<Collimator *>(coll.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
568

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

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

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

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

581
    double zstart = elptr->getZStart();
adelmann's avatar
adelmann committed
582
    *gmsg << "* Zstart= " << zstart << " [mm]" << endl;
583 584

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

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

590
    elptr->initialise(itsBunch);
gsell's avatar
gsell committed
591 592 593 594

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

596 597 598 599 600 601
    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 ;

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

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

657 658 659
void ParallelCyclotronTracker::visitOffset(const Offset & off) {
    if (opalRing_m == NULL)
        throw OpalException(
660 661
                            "ParallelCylcotronTracker::visitOffset",
                            "Attempt to place an offset when Ring not defined");
662 663
    Offset* offNonConst = const_cast<Offset*>(&off);
    offNonConst->updateGeometry(opalRing_m->getNextPosition(),
664
                                opalRing_m->getNextNormal());
665 666 667
    opalRing_m->appendElement(off);
}

gsell's avatar
gsell committed
668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689
/**
 *
 *
 * @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());
}

690

gsell's avatar
gsell committed
691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706
/**
 *
 *
 * @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) {
707
    *gmsg << "* -----------  Probe -------------------------------" << endl;
708 709
    Probe *elptr = dynamic_cast<Probe *>(prob.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
710

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

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

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

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

723
    double width = elptr->getWidth();
724
    *gmsg << "Width= " << width << " [mm]" << endl;
gsell's avatar
gsell committed
725 726 727


    // initialise, do nothing
728
    elptr->initialise(itsBunch);
gsell's avatar
gsell committed
729 730 731 732

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

734 735 736 737 738
    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
739 740

    // store probe parameters in the list
741
    buildupFieldList(BcParameter, ElementBase::PROBE, elptr);
gsell's avatar
gsell committed
742 743 744 745 746 747 748 749 750 751 752 753
}

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

754 755 756 757 758 759
void ParallelCyclotronTracker::visitSBend3D(const SBend3D &bend) {
    *gmsg << "Adding SBend3D" << endl;
    if (opalRing_m != NULL)
        opalRing_m->appendElement(bend);
    else
        throw OpalException("ParallelCyclotronTracker::visitSBend3D",
760
                            "Need to define a RINGDEFINITION to use SBend3D element");
761 762
}

763 764 765 766 767 768
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",
769
                            "Need to define a RINGDEFINITION to use VariableRFCavity element");
770 771
}

gsell's avatar
gsell committed
772 773 774 775 776 777 778 779
/**
 *
 *
 * @param as
 */
void ParallelCyclotronTracker::visitRFCavity(const RFCavity &as) {

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

781 782
    RFCavity *elptr = dynamic_cast<RFCavity *>(as.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
783 784 785 786 787 788 789

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

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

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

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

799
    std::string fmfn = elptr->getFieldMapFN();
gsell's avatar
gsell committed
800
    *gmsg << "* RF Field map file name= " << fmfn << endl;
801
    double angle = elptr->getAzimuth();
gsell's avatar
gsell committed
802 803
    *gmsg << "* Cavity azimuth position= " << angle << " [deg] " << endl;

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

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

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

813 814 815
    /*
      Setup time dependence and in case of no
      timedependence use a polynom with  a_0 = 1 and a_k = 0, k = 1,2,3.
816
    */
817

818 819 820 821
    std::shared_ptr<AbstractTimeDependence> freq_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> ampl_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> phase_atd = nullptr;

822 823 824 825 826
    std::vector<double>  unityVec;
    unityVec.push_back(1.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
827

828
    if (elptr->getFrequencyModelName() != "") {
829 830
        freq_atd = AbstractTimeDependence::getTimeDependence(elptr->getFrequencyModelName());
        *gmsg << "* Variable frequency RF Model name " << elptr->getFrequencyModelName() << endl;
831
    }
832 833
    else
        freq_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
834 835

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

    if (elptr->getPhaseModelName() != "") {
843 844
        phase_atd = AbstractTimeDependence::getTimeDependence(elptr->getPhaseModelName());
        *gmsg << "* Variable phase RF Model name " << elptr->getPhaseModelName() << endl;
845
    }
846 847
    else
        phase_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
848

gsell's avatar
gsell committed
849
    // read cavity voltage profile data from file.
850
    elptr->initialise(itsBunch, freq_atd, ampl_atd, phase_atd);
gsell's avatar
gsell committed
851 852 853 854

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

856 857
    BcParameter[0] = 0.001 * rmin;
    BcParameter[1] = 0.001 * rmax;
gsell's avatar
gsell committed
858 859 860
    BcParameter[2] = pdis;
    BcParameter[3] = angle;

861
    buildupFieldList(BcParameter, ElementBase::RFCAVITY, elptr);
gsell's avatar
gsell committed
862 863 864 865 866 867 868 869
}

/**
 *
 *
 * @param rfq
 */
void ParallelCyclotronTracker::visitRFQuadrupole(const RFQuadrupole &rfq) {
870
    *gmsg << "In RFQuadrupole; L = " << rfq.getElementLength() << " however the element is missing " << endl;
gsell's avatar
gsell committed
871 872 873 874 875 876 877 878 879
    myElements.push_back(dynamic_cast<RFQuadrupole *>(rfq.clone()));
}

/**
 *
 *
 * @param bend
 */
void ParallelCyclotronTracker::visitSBend(const SBend &bend) {
880
    *gmsg << "In SBend; L = " << bend.getElementLength() << " however the element is missing " << endl;
gsell's avatar
gsell committed
881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899
    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) {
900

901
    *gmsg << endl << "* -----------------------------  Septum ------------------------------- *" << endl;
gsell's avatar
gsell committed
902

903 904
    Septum *elptr = dynamic_cast<Septum *>(sept.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
905

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

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

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

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

918
    double width = elptr->getWidth();
919
    *gmsg << "Width = " << width << " [mm]" << endl;
gsell's avatar
gsell committed
920 921 922


    // initialise, do nothing
923
    elptr->initialise(itsBunch);
gsell's avatar
gsell committed
924 925 926 927

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

929 930 931