ParallelCyclotronTracker.h 16.8 KB
Newer Older
gsell's avatar
gsell committed
1 2
#ifndef OPAL_ParallelCyclotronTracker_HH
#define OPAL_ParallelCyclotronTracker_HH
3

gsell's avatar
gsell committed
4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
// ------------------------------------------------------------------------
// $RCSfile: ParallelCyclotronTracker.h,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.2.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: ParallelCyclotron
//
// ------------------------------------------------------------------------
//
// $Date: 2004/11/12 20:10:11 $
// $Author: adelmann $
//
// ------------------------------------------------------------------------

#include "Algorithms/Tracker.h"
#include "Structure/DataSink.h"
23
#include "AbsBeamline/ElementBase.h"
gsell's avatar
gsell committed
24 25 26 27 28
#include <vector>

class BMultipoleField;
class PartBunch;
class PlanarArcGeometry;
kraus's avatar
kraus committed
29
class Ring;
30
class SBend3D;
31 32
class VariableRFCavity;
class Offset;
gsell's avatar
gsell committed
33 34 35 36 37

// Class ParallelCyclotronTracker
// ------------------------------------------------------------------------
enum CyclOperationModeT {SINGLEP, MULTIP, TUNECALC};

38 39 40 41 42 43 44
struct CavityCrossData {
    RFCavity * cavity;
    double sinAzimuth;
    double cosAzimuth;
    double perpenDistance;
};

gsell's avatar
gsell committed
45 46 47 48
class ParallelCyclotronTracker: public Tracker {

public:

49
    typedef std::pair<double[8], Component *>      element_pair;
50
    typedef std::pair<ElementBase::ElementType, element_pair>        type_pair;
gsell's avatar
gsell committed
51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71
    typedef std::list<type_pair *>                 beamline_list;
    /// Constructor.
    //  The beam line to be tracked is "bl".
    //  The particle reference data are taken from "data".
    //  The particle bunch tracked is initially empty.
    //  If [b]revBeam[/b] is true, the beam runs from s = C to s = 0.
    //  If [b]revTrack[/b] is true, we track against the beam.
    explicit ParallelCyclotronTracker(const Beamline &bl, const PartData &data,
                                      bool revBeam, bool revTrack);

    /// Constructor.
    //  The beam line to be tracked is "bl".
    //  The particle reference data are taken from "data".
    //  The particle bunch tracked is taken from [b]bunch[/b].
    //  If [b]revBeam[/b] is true, the beam runs from s = C to s = 0.
    //  If [b]revTrack[/b] is true, we track against the beam.
    explicit ParallelCyclotronTracker(const Beamline &bl, PartBunch &bunch, DataSink &ds,
                                      const PartData &data, bool revBeam, bool revTrack, int maxSTEPS, int timeIntegrator);

    virtual ~ParallelCyclotronTracker();

kraus's avatar
kraus committed
72 73
    /// Apply the algorithm to an Ring
    virtual void visitRing(const Ring &ring);
74

gsell's avatar
gsell committed
75 76 77 78 79 80 81 82 83 84 85 86 87 88 89
    /// Apply the algorithm to a Cyclotorn
    virtual void visitCyclotron(const Cyclotron &cycl);

    /// Apply the algorithm to a RFCavity.
    virtual void visitRFCavity(const RFCavity &);

    /// Apply the algorithm to a BeamBeam.
    virtual void visitBeamBeam(const BeamBeam &);

    /// Apply the algorithm to a collimator.
    virtual void visitCollimator(const Collimator &);

    /// Apply the algorithm to a Corrector.
    virtual void visitCorrector(const Corrector &);

adelmann's avatar
adelmann committed
90 91 92
    /// Apply the algorithm to a Degrader
    virtual void visitDegrader(const Degrader &);

gsell's avatar
gsell committed
93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110
    /// Apply the algorithm to a Diagnostic.
    virtual void visitDiagnostic(const Diagnostic &);

    /// Apply the algorithm to a Drift.
    virtual void visitDrift(const Drift &);

    /// Apply the algorithm to a Lambertson.
    virtual void visitLambertson(const Lambertson &);

    /// Apply the algorithm to a Marker.
    virtual void visitMarker(const Marker &);

    /// Apply the algorithm to a Monitor.
    virtual void visitMonitor(const Monitor &);

    /// Apply the algorithm to a Multipole.
    virtual void visitMultipole(const Multipole &);

ext-rogers_c's avatar
ext-rogers_c committed
111 112 113
    /// Apply the algorithm to a MultipoleT
    virtual void visitMultipoleT (const MultipoleT &);
    
114 115 116
    /// Apply the algorithm to a Offset.
    virtual void visitOffset(const Offset &);

gsell's avatar
gsell committed
117 118 119 120 121 122 123 124 125 126 127 128
    /// Apply the algorithm to a Probe.
    virtual void visitProbe(const Probe &);

    /// Apply the algorithm to a RBend.
    virtual void visitRBend(const RBend &);

    /// Apply the algorithm to a RFQuadrupole.
    virtual void visitRFQuadrupole(const RFQuadrupole &);

    /// Apply the algorithm to a SBend.
    virtual void visitSBend(const SBend &);

129 130 131
    /// Apply the algorithm to a SBend3D.
    virtual void visitSBend3D(const SBend3D &);

gsell's avatar
gsell committed
132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149
    /// Apply the algorithm to a Separator.
    virtual void visitSeparator(const Separator &);

    /// Apply the algorithm to a Septum.
    virtual void visitSeptum(const Septum &);

    /// Apply the algorithm to a Solenoid.
    virtual void visitSolenoid(const Solenoid &);

    /// Apply the algorithm to a charge stripper.
    virtual void visitStripper(const Stripper &);

    /// Apply the algorithm to a ParallelPlate, it is empty for cyclotrontracker .
    virtual void visitParallelPlate(const ParallelPlate &);

    /// Apply the algorithm to a CyclotronValley.it is empty for cyclotrontracker .
    virtual void visitCyclotronValley(const CyclotronValley &);

150 151 152
    /// Apply the algorithm to a VariabelRFCavity.
    virtual void visitVariableRFCavity(const VariableRFCavity &cav);

gsell's avatar
gsell committed
153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174
    /// Apply the algorithm to the top-level beamline.
    //  overwrite the execute-methode from DefaultVisitor
    virtual void execute();

    /// Apply the algorithm to a beam line.
    //  overwrite the execute-methode from DefaultVisitor
    virtual void visitBeamline(const Beamline &);

    /// set total number of tracked bunches
    inline void setNumBunch(int n) { numBunch_m = n; }

    /// get total number of tracked bunches
    inline int  getNumBunch() { return numBunch_m; }

    /// set the working sub-mode for multi-bunch mode: "FORCE" or "AUTO"
    inline void  setMultiBunchMode(const int flag) {multiBunchMode_m = flag; }

    /// set last dumped step
    inline void setLastDumpedStep(const int para) {lastDumpedStep_m = para ; }

    /// set the control parameter for "AUTO" sub-mode
    inline void  setParaAutoMode(const double para) {CoeffDBunches_m = para; }
175

176
    inline void setPr(double x) {referencePr = x;}
177 178
    inline void setPt(double x) {referencePt = x;}
    inline void setPz(double x) {referencePz = x;}
179 180
    inline void setR(double x) {referenceR = x;}
    inline void setTheta(double x) {referenceTheta = x;}
181
    inline void setZ(double x) {referenceZ = x;}
182
    inline void setBeGa(double x) {bega = x;}
183 184 185
    inline void setPhi(double x) {referencePhi = x;}
    inline void setPsi(double x) {referencePsi = x;}
    inline void setPreviousH5Local(bool x) {previousH5Local = x;}
186 187 188

    void bgf_main_collision_test();
    void initializeBoundaryGeometry();
189

gsell's avatar
gsell committed
190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205
private:

    // Not implemented.
    ParallelCyclotronTracker();
    ParallelCyclotronTracker(const ParallelCyclotronTracker &);
    void operator=(const ParallelCyclotronTracker &);

    beamline_list FieldDimensions;
    std::list<Component *> myElements;
    int LastVisited;
    Beamline *itsBeamline;

    PartBunch *itsBunch;

    DataSink *itsDataSink;

206 207
    BoundaryGeometry *bgf_m;

gsell's avatar
gsell committed
208 209 210
    /// The maximal number of steps the system is integrated
    int maxSteps_m;

211
    /// The positive axes unit vectors
212 213 214
    static Vector_t const xaxis;
    static Vector_t const yaxis;
    static Vector_t const zaxis;
215

gsell's avatar
gsell committed
216 217
    /// The scale factor for dimensionless variables
    double scaleFactor_m;
Daniel Winklehner's avatar
Daniel Winklehner committed
218 219

    /// The reference variables
220
    double bega;
gsell's avatar
gsell committed
221 222
    double referenceR;
    double referenceTheta;
Daniel Winklehner's avatar
Daniel Winklehner committed
223
    double referenceZ = 0.0;
gsell's avatar
gsell committed
224 225 226

    double referencePr;
    double referencePt;
Daniel Winklehner's avatar
Daniel Winklehner committed
227
    double referencePz = 0.0;
gsell's avatar
gsell committed
228 229
    double referencePtot;

230 231 232
    double referencePsi;
    double referencePhi;

233 234
    bool spiral_flag = false;

235 236
    Vector_t PreviousMeanP;

237 238
    bool previousH5Local;

gsell's avatar
gsell committed
239 240
    double sinRefTheta_m;
    double cosRefTheta_m;
241

gsell's avatar
gsell committed
242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260
    /// The number of bunches specified in TURNS of RUN commond
    int numBunch_m;

    // 0 for single bunch (default),
    // 1 for FORCE,
    // 2 for AUTO
    int multiBunchMode_m;

    // control parameter for AUTO multi-bunch mode
    double CoeffDBunches_m;

    int lastDumpedStep_m;

    double PathLength_m;

    // the name of time integrator
    // The ID of time integrator
    // 0 --- RK-4(default)
    // 1 --- LF-2
261
    // 2 --- MTS
gsell's avatar
gsell committed
262 263 264 265
    int  timeIntegrator_m;

    void Tracker_LF();
    void Tracker_RK4();
266
    void Tracker_MTS();
gsell's avatar
gsell committed
267

268 269 270
    void Tracker_Generic();
    bool getFieldsAtPoint(const double &t, const size_t &Pindex, Vector_t &Efield, Vector_t &Bfield);

gsell's avatar
gsell committed
271
    /*
272
      Local Variables both used by the integration methods
gsell's avatar
gsell committed
273 274 275 276 277 278 279
    */

    Vector_t rold_m, pold_m, rnew_m, ptmp_m;

    long long step_m;
    long long restartStep0_m;

280 281 282 283
    int turnnumber_m;

    double const eta_m; // parameter for reset bin in multi-bunch run, todo: readin from inputfile

gsell's avatar
gsell committed
284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333
    // temporal 6 phase space varibles of particle [x,y,z,px,py,pz]. Unit: mm & dimensionless
    double variable_m[6];
    // temporal 3 real space varibles of particle ID=0 [x,y,z]. for tune with SC.  Unit: mm
    Vector_t variable_tune0_m;
    // temporal 3 real space varibles of particle ID=1 [x,y,z]. for tune with SC.  Unit: mm
    Vector_t variable_tune1_m;

    // vector of [angle, x, y] of SEO read in from external file for tune with SC. Unit : rad, mm
    std::vector<Vector_t> variable_SEO_m;

    // save initial phase space distribution (in global Cartesian frame ) for multi-bunch simultion. FixMe: not used
    Vector_t *initialR_m, *initialP_m;

    // record how many bunches has already been injected. ONLY FOR MPM
    int BunchCount_m;

    // decide how many energy bins. ONLY FOR MPM
    // For the time being, we set bin number equal to bunch number.
    int BinCount_m;

    // used for automatic injection in multi-bunch mode
    double RLastTurn_m , RThisTurn_m;

    // start time to record tune data
    double StartTime_m;

    // external field arrays for dumping
    Vector_t FDext_m[2], extE_m, extB_m;

    // mark the dumpstep to inject new bunch from here for AUTO mode of restart run of multibunch
    int backupDumpStep_m;

    const int myNode_m;
    const size_t initialLocalNum_m;
    const size_t initialTotalNum_m;

    std::ofstream outfTheta0_m;
    std::ofstream outfTheta1_m;
    std::ofstream outfTheta2_m;
    std::ofstream outfThetaEachTurn_m;


    //store the data of the beam which are required for injecting a new bunch for multibunch
    Ppos_t r_mb, p_mb;

    ParticleAttrib<double> q_mb;
    ParticleAttrib<double> m_mb;
    ParticleAttrib<short> ptype_mb;
    size_t npart_mb;

334
    void openFiles(std::string fn);
gsell's avatar
gsell committed
335 336 337 338 339 340 341 342
    void closeFiles();

    // Fringe fields for entrance and exit of magnetic elements.
    void applyEntranceFringe(double edge, double curve,
                             const BMultipoleField &field, double scale);
    void applyExitFringe(double edge, double curve,
                         const BMultipoleField &field, double scale);

343
    void buildupFieldList(double BcParameter[], ElementBase::ElementType elementType, Component *elptr);
gsell's avatar
gsell committed
344

345
    bool derivate(double *y, const double &t, double *yp, const size_t &Pindex);
gsell's avatar
gsell committed
346

347
    bool rk4(double x[], const double &t, const double &tau, const size_t &Pindex);
gsell's avatar
gsell committed
348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368

    // angle range [0~2PI) degree
    double calculateAngle(double x, double y);
    // angle range [-PI~PI) degree
    double calculateAngle2(double x, double y);

    bool readOneBunchFromFile(const size_t BeamCount);
    bool readOneBunch(const size_t BeamCount);
    void saveOneBunch();

    bool checkGapCross(Vector_t Rold, Vector_t Rnew, RFCavity * rfcavity, double &DistOld);
    bool RFkick(RFCavity * rfcavity, const double t, const double dt, const int Pindex);

    bool getTunes(std::vector<double> &t,  std::vector<double> &r,  std::vector<double> &z, int lastTurn, double &nur, double &nuz);

    IpplTimings::TimerRef IntegrationTimer_m;
    IpplTimings::TimerRef DumpTimer_m ;
    IpplTimings::TimerRef TransformTimer_m;
    IpplTimings::TimerRef BinRepartTimer_m;

    Vector_t calcMeanR() const;
369

gsell's avatar
gsell committed
370
    Vector_t calcMeanP() const;
371

gsell's avatar
gsell committed
372
    void repartition(); // Do repartition between nodes if step_m is multiple of Options::repartFreq
373 374

    // Transform the x- and y-parts of a particle attribute (position, momentum, fields) from the
gsell's avatar
gsell committed
375
    // global reference frame to the local reference frame.
376

gsell's avatar
gsell committed
377
    // phi is the angle of the bunch measured counter-clockwise from the positive x-axis.
378 379
    void globalToLocal(ParticleAttrib<Vector_t> & vectorArray,
                       double phi, Vector_t const translationToGlobal = 0);
380 381

    // Transform the x- and y-parts of a particle attribute (position, momentum, fields) from the
gsell's avatar
gsell committed
382
    // local reference frame to the global reference frame.
383
    void localToGlobal(ParticleAttrib<Vector_t> & vectorArray,
384
                       double phi, Vector_t const translationToGlobal = 0);
385

386
    // Overloaded version of globalToLocal using a quaternion for 3D rotation
387 388
    inline void globalToLocal(ParticleAttrib<Vector_t> & vectorArray,
                              Quaternion_t const quaternion,
389
                              Vector_t const meanR = Vector_t(0.0));
390

391
    // Overloaded version of localToGlobal using a quaternion for 3D rotation
392
    inline void localToGlobal(ParticleAttrib<Vector_t> & vectorArray,
393
                              Quaternion_t const quaternion,
394
                              Vector_t const meanR = Vector_t(0.0));
395

396
    // Overloaded version of globalToLocal using phi and theta for pseudo 3D rotation
397 398
    inline void globalToLocal(ParticleAttrib<Vector_t> & particleVectors,
                              double const phi, double const psi,
399
                              Vector_t const meanR = Vector_t(0.0));
400 401

    // Overloaded version of localToGlobal using phi and theta for pseudo 3D rotation
402 403
    inline void localToGlobal(ParticleAttrib<Vector_t> & particleVectors,
                              double const phi, double const psi,
404 405 406
                              Vector_t const meanR = Vector_t(0.0));

    // Overloaded version of globalToLocal using phi and theta for pseudo 3D rotation, single vector
407 408
    inline void globalToLocal(Vector_t & myVector,
                              double const phi, double const psi,
409 410 411
                              Vector_t const meanR = Vector_t(0.0));

    // Overloaded version of localToGlobal using phi and theta for pseudo 3D rotation, single vector
412 413
    inline void localToGlobal(Vector_t & myVector,
                              double const phi, double const psi,
414
                              Vector_t const meanR = Vector_t(0.0));
415 416

    // Rotate the particles by an angle and axis defined in the quaternion (4-vector w,x,y,z)
adelmann's avatar
adelmann committed
417
    inline void rotateWithQuaternion(ParticleAttrib<Vector_t> & vectorArray, Quaternion_t const quaternion);
418 419

    // Returns the quaternion of the rotation from vector u to vector v
adelmann's avatar
adelmann committed
420
    inline void getQuaternionTwoVectors(Vector_t u, Vector_t v, Quaternion_t & quaternion);
421 422

    // Normalization of a quaternion
adelmann's avatar
adelmann committed
423
    inline void normalizeQuaternion(Quaternion_t & quaternion);
424

425
    // Normalization of a quaternion
426
    inline void normalizeVector(Vector_t & vector);
427 428 429 430

    // Some rotations
    inline void rotateAroundZ(ParticleAttrib<Vector_t> & particleVectors, double const phi);
    inline void rotateAroundX(ParticleAttrib<Vector_t> & particleVectors, double const psi);
431 432
    inline void rotateAroundZ(Vector_t & myVector, double const phi);
    inline void rotateAroundX(Vector_t & myVector, double const psi);
433

434
    // Push particles for time h.
gsell's avatar
gsell committed
435 436 437 438 439
    // Apply effects of RF Gap Crossings.
    // Update time and path length.
    // Unit assumptions: [itsBunch->R] = m, [itsBunch->P] = 1, [h] = s, [c] = m/s, [itsBunch->getT()] = s
    void push(double h);

440 441 442 443 444 445 446
    // Kick particles for time h
    // The fields itsBunch->Bf, itsBunch->Ef are used to calculate the forces
    void kick(double h);

    // Apply the trilogy half push - kick - half push,
    // considering only external fields
    void borisExternalFields(double h);
447 448

    // apply the plugin elements: probe, collimator, stripper, septum
449
    void applyPluginElements(const double dt);
450

451 452
    // destroy particles if they are marked as Bin=-1 in the plugin elements or out of global apeture
    bool deleteParticle();
453

454 455 456 457 458 459
    std::ofstream outfTrackOrbit_m;

    void initTrackOrbitFile();

    void singleParticleDump();

460 461 462
    void bunchDumpStatData();

    void bunchDumpPhaseSpaceData();
463

464 465
    void evaluateSpaceChargeField();

466 467
    void initDistInGlobalFrame();

Matthias Toggweiler's avatar
Matthias Toggweiler committed
468 469
    void checkNumPart(std::string s);

kraus's avatar
kraus committed
470 471
    // we store a pointer explicitly to the Ring
    Ring* opalRing_m;
472

kraus's avatar
kraus committed
473
    // If Ring is defined take the harmonic number from Ring; else use
474 475 476
    // cyclotron
    double getHarmonicNumber() const;

gsell's avatar
gsell committed
477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503
};

/**
 *
 *
 * @param x
 * @param y
 *
 * @return angle range [0~2PI) degree
 */
inline
double ParallelCyclotronTracker::calculateAngle(double x, double y) {
    double thetaXY = atan2(y, x);

    if (thetaXY < 0) return thetaXY + Physics::two_pi;
    return thetaXY;
}

/**
 *
 *
 * @param x
 * @param y
 *
 * @return angle range [-PI~PI) degree
 */
inline
504
double ParallelCyclotronTracker::calculateAngle2(double x, double y)
gsell's avatar
gsell committed
505 506
{ return atan2(y,x); }

winklehner_d's avatar
winklehner_d committed
507
#endif // OPAL_ParallelCyclotronTracker_HH