ParallelCyclotronTracker.h 18.5 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
// ------------------------------------------------------------------------
// $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"
22
#include "AbsBeamline/ElementBase.h"
gsell's avatar
gsell committed
23
#include <vector>
24
#include <tuple>
gsell's avatar
gsell committed
25

26
#include "Steppers/Steppers.h"
gsell's avatar
gsell committed
27

frey_m's avatar
frey_m committed
28 29
#include "Structure/MultiBunchDump.h"

30 31
class DataSink;

frey_m's avatar
frey_m committed
32 33
template <class T, unsigned Dim>
class PartBunchBase;
gsell's avatar
gsell committed
34 35 36 37

// Class ParallelCyclotronTracker
// ------------------------------------------------------------------------

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
class ParallelCyclotronTracker: public Tracker {

public:
48

49
    enum class MODE {
50 51 52 53 54
        UNDEFINED = -1,
        SINGLE = 0,
        SEO = 1,
        BUNCH = 2
    };
55

56 57 58 59 60 61
    // multi-bunch modes
    enum class MB_MODE {
        NONE   = 0,
        FORCE  = 1,
        AUTO   = 2
    };
62

frey_m's avatar
frey_m committed
63 64 65 66 67
    // multi-bunch binning type
    enum class MB_BINNING {
        GAMMA = 0,
        BUNCH = 1
    };
gsell's avatar
gsell committed
68

69 70
    typedef std::vector<double> dvector_t;
    typedef std::vector<int> ivector_t;
71
    typedef std::pair<double[8], Component *>      element_pair;
72
    typedef std::pair<ElementBase::ElementType, element_pair>        type_pair;
gsell's avatar
gsell committed
73 74 75 76 77 78 79 80
    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 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.
81 82 83 84
    ParallelCyclotronTracker(const Beamline &bl, PartBunchBase<double, 3> *bunch, DataSink &ds,
                             const PartData &data, bool revBeam,
                             bool revTrack, int maxSTEPS,
                             int timeIntegrator, int numBunch);
gsell's avatar
gsell committed
85 86 87

    virtual ~ParallelCyclotronTracker();

kraus's avatar
kraus committed
88 89
    /// Apply the algorithm to an Ring
    virtual void visitRing(const Ring &ring);
90

gsell's avatar
gsell committed
91 92 93 94 95 96 97 98 99 100
    /// 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.
101
    virtual void visitCCollimator(const CCollimator &);
gsell's avatar
gsell committed
102 103 104 105

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

adelmann's avatar
adelmann committed
106 107 108
    /// Apply the algorithm to a Degrader
    virtual void visitDegrader(const Degrader &);

gsell's avatar
gsell committed
109 110 111 112 113 114
    /// Apply the algorithm to a Diagnostic.
    virtual void visitDiagnostic(const Diagnostic &);

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

115 116 117
    /// Apply the algorithm to a flexible collimator
    virtual void visitFlexibleCollimator(const FlexibleCollimator &);

gsell's avatar
gsell committed
118 119 120 121 122 123 124 125 126 127 128 129
    /// 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
130 131
    /// Apply the algorithm to a MultipoleT
    virtual void visitMultipoleT (const MultipoleT &);
132 133 134 135 136 137 138 139 140

    /// Apply the algorithm to a MultipoleTStraight
    virtual void visitMultipoleTStraight (const MultipoleTStraight &);

    /// Apply the algorithm to a MultipoleTCurvedConstRadius
    virtual void visitMultipoleTCurvedConstRadius (const MultipoleTCurvedConstRadius &);

    /// Apply the algorithm to a MultipoleTCurvedVarRadius
    virtual void visitMultipoleTCurvedVarRadius (const MultipoleTCurvedVarRadius &);
141

142 143 144
    /// Apply the algorithm to a Offset.
    virtual void visitOffset(const Offset &);

gsell's avatar
gsell committed
145 146 147 148 149 150 151 152 153 154 155 156
    /// 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 &);

157 158 159
    /// Apply the algorithm to a SBend3D.
    virtual void visitSBend3D(const SBend3D &);

ext-rogers_c's avatar
ext-rogers_c committed
160 161
    /// Apply the algorithm to a ScalingFFAMagnet.
    virtual void visitScalingFFAMagnet(const ScalingFFAMagnet &bend);
162

gsell's avatar
gsell committed
163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
    /// 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 &);

181 182 183
    /// Apply the algorithm to a VariabelRFCavity.
    virtual void visitVariableRFCavity(const VariableRFCavity &cav);

184 185 186 187
    /// Apply the algorithm to a VariabelRFCavity.
    virtual void visitVariableRFCavityFringeField
                                      (const VariableRFCavityFringeField &cav);

gsell's avatar
gsell committed
188 189 190 191 192 193 194 195 196 197 198 199 200 201 202
    /// 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"
203
    void setMultiBunchMode(const std::string& mbmode);
204

frey_m's avatar
frey_m committed
205 206
    // set binning type
    void setMultiBunchBinning(std::string binning);
207

frey_m's avatar
frey_m committed
208 209
    /// set the scale for binning in multi-bunch mode
    void setMultiBunchEta(const double& eta) { eta_m = eta; };
gsell's avatar
gsell committed
210 211 212 213 214

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

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

217
    ///@{ Method for restart
218
    inline void setPr(double x) {referencePr = x;}
219 220
    inline void setPt(double x) {referencePt = x;}
    inline void setPz(double x) {referencePz = x;}
221 222
    inline void setR(double x) {referenceR = x;}
    inline void setTheta(double x) {referenceTheta = x;}
223
    inline void setZ(double x) {referenceZ = x;}
224
    inline void setBeGa(double x) {bega = x;}
225 226 227
    inline void setPhi(double x) {referencePhi = x;}
    inline void setPsi(double x) {referencePsi = x;}
    inline void setPreviousH5Local(bool x) {previousH5Local = x;}
228
    ///@}
229 230
    void bgf_main_collision_test();
    void initializeBoundaryGeometry();
231

gsell's avatar
gsell committed
232 233 234 235 236 237 238 239 240 241 242 243
private:

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

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

    DataSink *itsDataSink;
244

frey_m's avatar
frey_m committed
245
    std::unique_ptr<MultiBunchDump> itsMBDump_m;
gsell's avatar
gsell committed
246

247 248
    BoundaryGeometry *bgf_m;

gsell's avatar
gsell committed
249 250 251
    /// The maximal number of steps the system is integrated
    int maxSteps_m;

252
    /// The positive axes unit vectors
253 254 255
    static Vector_t const xaxis;
    static Vector_t const yaxis;
    static Vector_t const zaxis;
256

Daniel Winklehner's avatar
Daniel Winklehner committed
257
    /// The reference variables
258
    double bega;
gsell's avatar
gsell committed
259 260
    double referenceR;
    double referenceTheta;
Daniel Winklehner's avatar
Daniel Winklehner committed
261
    double referenceZ = 0.0;
gsell's avatar
gsell committed
262 263 264

    double referencePr;
    double referencePt;
Daniel Winklehner's avatar
Daniel Winklehner committed
265
    double referencePz = 0.0;
gsell's avatar
gsell committed
266 267
    double referencePtot;

268 269 270
    double referencePsi;
    double referencePhi;

271 272
    bool spiral_flag = false;

273 274
    Vector_t PreviousMeanP;

275 276
    bool previousH5Local;

gsell's avatar
gsell committed
277 278
    double sinRefTheta_m;
    double cosRefTheta_m;
279

gsell's avatar
gsell committed
280 281 282 283 284 285
    /// The number of bunches specified in TURNS of RUN commond
    int numBunch_m;

    // 0 for single bunch (default),
    // 1 for FORCE,
    // 2 for AUTO
286
    MB_MODE multiBunchMode_m;
287

frey_m's avatar
frey_m committed
288 289 290
    // 0 for GAMMA (default),
    // 1 for BUNCH
    MB_BINNING binningType_m;
gsell's avatar
gsell committed
291 292 293 294 295 296 297 298

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

    int lastDumpedStep_m;

    double PathLength_m;

299
    void MtsTracker();
gsell's avatar
gsell committed
300

301
    void GenericTracker();
302 303
    bool getFieldsAtPoint(const double &t, const size_t &Pindex, Vector_t &Efield, Vector_t &Bfield);

gsell's avatar
gsell committed
304
    /*
305
      Local Variables both used by the integration methods
gsell's avatar
gsell committed
306 307 308 309 310
    */

    long long step_m;
    long long restartStep0_m;

311 312
    int turnnumber_m;

frey_m's avatar
frey_m committed
313
    double eta_m; // parameter for reset bin in multi-bunch run
314

gsell's avatar
gsell committed
315 316 317 318 319 320 321 322 323 324 325 326 327
    // record how many bunches has already been injected. ONLY FOR MPM
    int BunchCount_m;

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

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

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

328 329 330 331 332 333 334 335 336 337
    /// output coordinates at different azimuthal angles and one after every turn
    std::vector<std::ofstream> outfTheta_m;
    /// the different azimuthal angles for the outfTheta_m output files
    std::vector<double> azimuth_angle_m;
    ///@ open / close output coordinate files
    void openFiles(std::string fn);
    void closeFiles();
    ///@}
    /// output file for six dimensional phase space
    std::ofstream outfTrackOrbit_m;
gsell's avatar
gsell committed
338 339

    //store the data of the beam which are required for injecting a new bunch for multibunch
340 341
    /// filename
    std::string onebunch_m;
342

gsell's avatar
gsell committed
343

344
    void buildupFieldList(double BcParameter[], ElementBase::ElementType elementType, Component *elptr);
gsell's avatar
gsell committed
345 346 347 348 349 350 351 352

    // 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);
    void saveOneBunch();
353

frey_m's avatar
frey_m committed
354
    void updateParticleBins_m();
gsell's avatar
gsell committed
355 356 357 358

    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);

359
    bool getTunes(dvector_t &t,  dvector_t &r,  dvector_t &z, int lastTurn, double &nur, double &nuz);
gsell's avatar
gsell committed
360 361 362 363 364 365 366

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

    Vector_t calcMeanR() const;
367

gsell's avatar
gsell committed
368
    Vector_t calcMeanP() const;
369

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

438
    // Kick particles for time h
frey_m's avatar
frey_m committed
439
    // The fields itsBunch_m->Bf, itsBunch_m->Ef are used to calculate the forces
440 441 442 443 444
    void kick(double h);

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

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

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

452 453 454 455
    void initTrackOrbitFile();

    void singleParticleDump();

456
    void bunchDumpStatData();
457

458 459
    // @param azimuth (global) [deg] of dump
    void bunchDumpStatDataPerBin(const double& azimuth);
460 461

    void bunchDumpPhaseSpaceData();
462

463 464
    void evaluateSpaceChargeField();

465 466
    void initDistInGlobalFrame();

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

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

472

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

477 478 479 480
    typedef std::function<bool(const double&,
                               const size_t&,
                               Vector_t&,
                               Vector_t&)> function_t;
481

482
    std::unique_ptr< Stepper<function_t> > itsStepper_mp;
483

484 485 486 487 488 489
    struct settings {
        int scSolveFreq;
        int stepsPerTurn;
        double deltaTheta;
        int stepsNextCheck;
    } setup_m;
490

491
    MODE mode_m;
492

493
    stepper::INTEGRATOR stepper_m;
494

495
    void update_m(double& t, const double& dt, const bool& dumpEachTurn);
496

497 498 499
    /*!
     * @returns the time t [ns], time step dt [ns] and the azimuth angle [rad]
     */
500
    std::tuple<double, double, double> initializeTracking_m();
501

502 503 504 505
    void finalizeTracking_m(dvector_t& Ttime,
                            dvector_t& Tdeltr,
                            dvector_t& Tdeltz,
                            ivector_t& TturnNumber);
506

507
    void seoMode_m(double& t, const double dt, bool& dumpEachTurn,
508 509
                   dvector_t& Ttime, dvector_t& Tdeltr,
                   dvector_t& Tdeltz, ivector_t& TturnNumber);
510

511
    void singleMode_m(double& t, const double dt, bool& dumpEachTurn, double& oldReferenceTheta);
512

513
    void bunchMode_m(double& t, const double dt, bool& dumpEachTurn);
514

515 516
    void gapCrossKick_m(size_t i, double t, double dt,
                        const Vector_t& Rold, const Vector_t& Pold);
517 518


519 520 521 522 523
    inline void dumpAzimuthAngles_m(const double& t,
                                    const Vector_t& R,
                                    const Vector_t& P,
                                    const double& oldReferenceTheta,
                                    const double& temp_meanTheta);
524

525 526 527 528 529
    inline void dumpThetaEachTurn_m(const double& t,
                                    const Vector_t& R,
                                    const Vector_t& P,
                                    const double& temp_meanTheta,
                                    bool& dumpEachTurn);
530

531
    void computeSpaceChargeFields_m();
532

533 534 535 536
    bool computeExternalFields_m(const size_t& i,
                                 const double& t,
                                 Vector_t& Efield,
                                 Vector_t& Bfield);
537

538
    void injectBunch_m(bool& flagTransition);
539

gsell's avatar
gsell committed
540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562
};

/**
 * @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
563
double ParallelCyclotronTracker::calculateAngle2(double x, double y)
gsell's avatar
gsell committed
564 565
{ return atan2(y,x); }

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