ParallelCyclotronTracker.cpp 155 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-9;  // 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
    // This has to match the dT in the rk4 pusher
    double dtime = itsBunch->getdT() * getHarmonicNumber();
227

228 229
    int triId = 0;
    for(size_t i = 0; i < itsBunch->getLocalNum(); i++) {
230 231
        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);
232 233 234
        if(res >= 0) {
            itsBunch->Bin[i] = -1;
        }
235
    }
236 237 238
}


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

246
    std::string  SfileName2 = SfileName + std::string("-Angle0.dat");
gsell's avatar
gsell committed
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.) For high currentst 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
        // Get reference values from cyclotron element 
        // For now, these are still stored in mm. should be the only ones. -DW
388 389 390 391
        referenceR     = elptr->getRinit();
        referenceTheta = elptr->getPHIinit();
        referenceZ     = elptr->getZinit();
        referencePr    = elptr->getPRinit();
392 393
        referencePz    = elptr->getPZinit();

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

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

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

        if(insqrt < 0) {

            if(insqrt > -1.0e-10) {

	        referencePt = 0.0;

            } else {

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

        } else {

            referencePt = sqrt(insqrt);
        }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

525
    // Read in cyclotron field maps (midplane + 3D fields if desired).
526
    elptr->initialise(itsBunch, fieldflag, elptr->getBScale());
gsell's avatar
gsell committed
527 528

    double BcParameter[8];
529

530
    for(int i = 0; i < 8; i++)
531
        BcParameter[i] = 0.0;
532

533 534
    BcParameter[0] = 0.001 * elptr->getRmin();
    BcParameter[1] = 0.001 * elptr->getRmax();
gsell's avatar
gsell committed
535

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

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

555
    *gmsg << "* --------- Collimator -----------------------------" << endl;
gsell's avatar
gsell committed
556

557 558
    Collimator* elptr = dynamic_cast<Collimator *>(coll.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
559

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

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

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

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

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

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

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

581
    elptr->initialise(itsBunch);
gsell's avatar
gsell committed
582 583 584 585

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

587 588 589 590 591 592
    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 ;

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

/**
 *
 *
 * @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
606 607 608 609 610 611 612 613 614 615 616 617
/**
 *
 *
 * @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
618 619 620 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
/**
 *
 *
 * @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()));
}

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

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

681

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

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

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

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

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

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


    // initialise, do nothing
719
    elptr->initialise(itsBunch);
gsell's avatar
gsell committed
720 721 722 723

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

725 726 727 728 729
    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
730 731

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

810 811 812 813
    std::shared_ptr<AbstractTimeDependence> freq_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> ampl_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> phase_atd = nullptr;

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

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

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

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

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

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

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

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

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

/**
 *
 *
 * @param bend
 */
void ParallelCyclotronTracker::visitSBend(const SBend &bend) {
872
    *gmsg << "In SBend; L = " << bend.getElementLength() << " however the element is missing " << endl;
gsell's avatar
gsell committed
873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891
    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) {
892

893
    *gmsg << endl << "* -----------------------------  Septum ------------------------------- *" << endl;
gsell's avatar
gsell committed
894

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

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

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

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

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

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


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

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

921 922 923 924 925
    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
926 927

    // store septum parameters in the list
928
    buildupFieldList(BcParameter, ElementBase::SEPTUM, elptr);
gsell's avatar
gsell committed
929 930 931