ParallelSliceTracker.h 11.2 KB
Newer Older
gsell's avatar
gsell committed
1 2 3 4 5 6
#ifndef OPAL_ParallelSliceTracker_HH
#define OPAL_ParallelSliceTracker_HH

#include <vector>
#include <list>

7 8
class EnvelopeBunch;

9
#include "Algorithms/Vektor.h"
gsell's avatar
gsell committed
10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55
#include "Algorithms/Tracker.h"
#include "Structure/DataSink.h"
#include "Utilities/Options.h"

#include "Physics/Physics.h"

#include "AbsBeamline/AlignWrapper.h"
#include "AbsBeamline/BeamBeam.h"
#include "AbsBeamline/Collimator.h"
#include "AbsBeamline/Corrector.h"
#include "AbsBeamline/Diagnostic.h"
#include "AbsBeamline/Drift.h"
#include "AbsBeamline/ElementBase.h"
#include "AbsBeamline/Lambertson.h"
#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/TravelingWave.h"
#include "AbsBeamline/RFQuadrupole.h"
#include "AbsBeamline/SBend.h"
#include "AbsBeamline/Separator.h"
#include "AbsBeamline/Septum.h"
#include "AbsBeamline/Solenoid.h"
#include "AbsBeamline/ParallelPlate.h"
#include "AbsBeamline/CyclotronValley.h"

#include "Algorithms/ParallelTTracker.h"

#include "Beamlines/Beamline.h"
#include "Elements/OpalBeamline.h"

class ParallelSliceTracker: public Tracker {

public:

    //  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 ParallelSliceTracker(const Beamline &bl, const PartData &data,
                                  bool revBeam, bool revTrack);

    //  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.
56 57 58 59
    explicit ParallelSliceTracker(const Beamline &bl, EnvelopeBunch &bunch,
                                  DataSink &ds, const PartData &data,
                                  bool revBeam, bool revTrack, int maxSTEPS,
                                  double zstop, ParallelTTracker &mySlApTracker);
gsell's avatar
gsell committed
60 61 62 63 64 65 66

    virtual ~ParallelSliceTracker();

    virtual void visitAlignWrapper(const AlignWrapper &);
    virtual void visitBeamBeam(const BeamBeam &);
    virtual void visitCollimator(const Collimator &);
    virtual void visitCorrector(const Corrector &);
adelmann's avatar
adelmann committed
67
    virtual void visitDegrader(const Degrader &);
gsell's avatar
gsell committed
68 69 70 71 72
    virtual void visitDiagnostic(const Diagnostic &);
    virtual void visitDrift(const Drift &);
    virtual void visitLambertson(const Lambertson &);
    virtual void visitMarker(const Marker &);
    virtual void visitMonitor(const Monitor &);
73
    virtual void visitMultipole(const Multipole &);
gsell's avatar
gsell committed
74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91
    virtual void visitProbe(const Probe &);
    virtual void visitRBend(const RBend &);
    virtual void visitRFCavity(const RFCavity &);
    virtual void visitTravelingWave(const TravelingWave &);
    virtual void visitRFQuadrupole(const RFQuadrupole &);
    virtual void visitSBend(const SBend &);
    virtual void visitSeparator(const Separator &);
    virtual void visitSeptum(const Septum &);
    virtual void visitSolenoid(const Solenoid &);
    virtual void visitParallelPlate(const ParallelPlate &);
    virtual void visitCyclotronValley(const CyclotronValley &);

    virtual void execute();

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

92

gsell's avatar
gsell committed
93 94 95 96 97 98 99 100 101 102
private:

    ParallelSliceTracker();
    ParallelSliceTracker(const ParallelSliceTracker &);
    void operator=(const ParallelSliceTracker &);

    void updateRFElement(string elName, double maxPhi);
    void updateAllRFElements();
    double getCavityPhase(FieldList cav, string name);

103 104
    double currentSimulationTime_m;
    bool globalEOL_m;
gsell's avatar
gsell committed
105

106
    std::unique_ptr<OpalBeamline> itsOpalBeamline_m;
gsell's avatar
gsell committed
107

108
    EnvelopeBunch *itsBunch_m;
gsell's avatar
gsell committed
109 110 111

    ParallelTTracker *mySlApTracker_m;

112
    DataSink *itsDataSink_m;
gsell's avatar
gsell committed
113 114

    /// The maximal number of steps the system is integrated
115
    unsigned long long maxSteps_m;
gsell's avatar
gsell committed
116 117 118 119 120 121 122 123 124

    double zstop_m;

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

125 126
    void FieldsOutput(double z, double Ex, double Ey, double Ez,
                      double Bx, double By, double Bz);
gsell's avatar
gsell committed
127 128 129 130 131

    void kickParticles();

    void updateReferenceParticle();
    void updateSpaceOrientation(const bool &move = false);
132 133
    void kickReferenceParticle(const Vector_t &externalE,
                               const Vector_t &externalB);
gsell's avatar
gsell committed
134 135 136
    void writePhaseSpace(const long long step, const double &sposRef);
    void writeLastStepPhaseSpace(const long long step, const double &sposRef);

137 138 139 140 141 142 143 144 145 146 147 148 149
    void prepareSections();
    void doAutoPhasing();
    void timeIntegration();
    void computeExternalFields();
    void switchElements(double scaleMargin = 3.0);
    void computeSpaceChargeFields();
    void dumpStats(long long step);
    bool hasEndOfLineReached();
    void handleRestartRun();
    void setTime();
    void setLastStep();
    void dumpPhaseSpaceOnScan();

gsell's avatar
gsell committed
150 151 152 153 154 155 156 157 158 159 160 161
    IpplTimings::TimerRef timeIntegrationTimer1_m;
    IpplTimings::TimerRef timeIntegrationTimer2_m;
    IpplTimings::TimerRef timeFieldEvaluation_m ;
    IpplTimings::TimerRef BinRepartTimer_m;
    IpplTimings::TimerRef WakeFieldTimer_m;
};

inline void ParallelSliceTracker::visitBeamline(const Beamline &bl) {
    itsBeamline_m.iterate(*dynamic_cast<BeamlineVisitor *>(this), false);
}

inline void ParallelSliceTracker::visitAlignWrapper(const AlignWrapper &wrap) {
162
    itsOpalBeamline_m->visit(wrap, *this, itsBunch_m);
gsell's avatar
gsell committed
163 164 165
}

inline void ParallelSliceTracker::visitBeamBeam(const BeamBeam &bb) {
166
    itsOpalBeamline_m->visit(bb, *this, itsBunch_m);
gsell's avatar
gsell committed
167 168 169
}

inline void ParallelSliceTracker::visitCollimator(const Collimator &coll) {
170
    itsOpalBeamline_m->visit(coll, *this, itsBunch_m);
gsell's avatar
gsell committed
171 172 173 174
}


inline void ParallelSliceTracker::visitCorrector(const Corrector &corr) {
175
    itsOpalBeamline_m->visit(corr, *this, itsBunch_m);
gsell's avatar
gsell committed
176 177
}

adelmann's avatar
adelmann committed
178 179 180
inline void ParallelSliceTracker::visitDegrader(const Degrader &deg) {
    itsOpalBeamline_m->visit(deg, *this, itsBunch_m);
}
gsell's avatar
gsell committed
181 182

inline void ParallelSliceTracker::visitDiagnostic(const Diagnostic &diag) {
183
    itsOpalBeamline_m->visit(diag, *this, itsBunch_m);
gsell's avatar
gsell committed
184 185 186 187
}


inline void ParallelSliceTracker::visitDrift(const Drift &drift) {
188
    itsOpalBeamline_m->visit(drift, *this, itsBunch_m);
gsell's avatar
gsell committed
189 190 191 192
}


inline void ParallelSliceTracker::visitLambertson(const Lambertson &lamb) {
193
    itsOpalBeamline_m->visit(lamb, *this, itsBunch_m);
gsell's avatar
gsell committed
194 195 196 197
}


inline void ParallelSliceTracker::visitMarker(const Marker &marker) {
198
    itsOpalBeamline_m->visit(marker, *this, itsBunch_m);
gsell's avatar
gsell committed
199 200 201 202
}


inline void ParallelSliceTracker::visitMonitor(const Monitor &mon) {
203
    itsOpalBeamline_m->visit(mon, *this, itsBunch_m);
gsell's avatar
gsell committed
204 205 206 207
}


inline void ParallelSliceTracker::visitMultipole(const Multipole &mult) {
208
    itsOpalBeamline_m->visit(mult, *this, itsBunch_m);
gsell's avatar
gsell committed
209 210 211 212
}


inline void ParallelSliceTracker::visitProbe(const Probe &prob) {
213
    itsOpalBeamline_m->visit(prob, *this, itsBunch_m);
gsell's avatar
gsell committed
214 215 216 217
}


inline void ParallelSliceTracker::visitRBend(const RBend &bend) {
218
    itsOpalBeamline_m->visit(bend, *this, itsBunch_m);
gsell's avatar
gsell committed
219 220 221 222
}


inline void ParallelSliceTracker::visitRFCavity(const RFCavity &as) {
223
    itsOpalBeamline_m->visit(as, *this, itsBunch_m);
gsell's avatar
gsell committed
224 225 226
}

inline void ParallelSliceTracker::visitTravelingWave(const TravelingWave &as) {
227
    itsOpalBeamline_m->visit(as, *this, itsBunch_m);
gsell's avatar
gsell committed
228 229 230 231
}


inline void ParallelSliceTracker::visitRFQuadrupole(const RFQuadrupole &rfq) {
232
    itsOpalBeamline_m->visit(rfq, *this, itsBunch_m);
gsell's avatar
gsell committed
233 234 235
}

inline void ParallelSliceTracker::visitSBend(const SBend &bend) {
236
    itsOpalBeamline_m->visit(bend, *this, itsBunch_m);
gsell's avatar
gsell committed
237 238 239 240
}


inline void ParallelSliceTracker::visitSeparator(const Separator &sep) {
241
    itsOpalBeamline_m->visit(sep, *this, itsBunch_m);
gsell's avatar
gsell committed
242 243 244 245
}


inline void ParallelSliceTracker::visitSeptum(const Septum &sept) {
246
    itsOpalBeamline_m->visit(sept, *this, itsBunch_m);
gsell's avatar
gsell committed
247 248 249 250
}


inline void ParallelSliceTracker::visitSolenoid(const Solenoid &solenoid) {
251
    itsOpalBeamline_m->visit(solenoid, *this, itsBunch_m);
gsell's avatar
gsell committed
252 253 254 255 256 257 258 259 260 261 262 263 264 265
}

inline void ParallelSliceTracker::visitParallelPlate(const ParallelPlate &pplate) {
    //do nothing.
}

inline void ParallelSliceTracker::visitCyclotronValley(const CyclotronValley &cv) {
    // Do nothing here.
}

inline void ParallelSliceTracker::kickParticles() {
}

inline void ParallelSliceTracker::updateReferenceParticle() {
266
    itsBunch_m->calcBeamParameters();
gsell's avatar
gsell committed
267 268 269 270

}

inline void ParallelSliceTracker::updateSpaceOrientation(const bool &move) {
271
    itsBunch_m->calcBeamParameters();
gsell's avatar
gsell committed
272 273 274 275 276 277 278 279
}

inline void ParallelSliceTracker::kickReferenceParticle(const Vector_t &externalE, const Vector_t &externalB) {
}

inline void ParallelSliceTracker::writePhaseSpace(const long long step, const double &sposRef) {
    Inform msg("ParallelSliceTracker");
    Vector_t externalE, externalB;
280
    Vector_t FDext[6];  // = {BHead, EHead, BRef, ERef, BTail, ETail}.
gsell's avatar
gsell committed
281

282 283 284 285 286
    /*
     * Sample fields at (xmin, ymin, zmin), (xmax, ymax, zmax) and the
     * centroid location. We are sampling the electric and magnetic fields at
     * the back, front and center of the beam.
     */
gsell's avatar
gsell committed
287
    Vector_t rmin, rmax;
288
    itsBunch_m->get_bounds(rmin, rmax);
gsell's avatar
gsell committed
289 290 291 292 293 294 295 296 297

    Vector_t pos[3];
    pos[0] = Vector_t(rmax(0), rmax(1), rmax(2));
    pos[1] = Vector_t(0.0, 0.0, sposRef);
    pos[2] = Vector_t(rmin(0), rmin(1), rmin(2));

    for(int k = 0; k < 3; ++k) {
        externalB = Vector_t(0.0);
        externalE = Vector_t(0.0);
298 299 300
        double bunch_time = itsBunch_m->getT() - 0.5 * itsBunch_m->getdT();
        itsOpalBeamline_m->getFieldAt(pos[k], itsBunch_m->get_rmean(),
                                      bunch_time, externalE, externalB);
gsell's avatar
gsell committed
301 302 303 304 305 306

        FDext[2*k]   = externalB;
        FDext[2*k+1] = externalE * 1e-6;
    }

    if(step % Options::psDumpFreq == 0) {
307 308 309 310 311
        //itsDataSink->stashPhaseSpaceEnvelope(*itsBunch_m, FDext, rmax(2), sposRef, rmin(2));
        itsDataSink_m->writePhaseSpaceEnvelope(*itsBunch_m, FDext,
                                             rmax(2), sposRef, rmin(2));
        INFOMSG("* Wrote beam phase space." << endl);
        msg     << *itsBunch_m << endl;
gsell's avatar
gsell committed
312 313 314
    }

    if(step % Options::statDumpFreq == 0) {
315 316
        itsDataSink_m->writeStatData(*itsBunch_m, FDext,
                                   rmax(2), sposRef, rmin(2));
gsell's avatar
gsell committed
317 318 319 320 321 322
        INFOMSG("* Wrote beam statistics." << endl);
    }
}

inline void ParallelSliceTracker::writeLastStepPhaseSpace(const long long step, const double &sposRef) {
    Inform msg("ParallelSliceTracker");
323
    if(itsBunch_m->isValid_m) {
gsell's avatar
gsell committed
324 325 326
        Vector_t externalE, externalB;
        Vector_t FDext[6];
        Vector_t rmin, rmax;
327
        itsBunch_m->get_bounds(rmin, rmax);
gsell's avatar
gsell committed
328 329 330 331 332 333 334 335 336 337 338 339 340

        Vector_t pos[3];
        pos[0] = Vector_t(rmax(0), rmax(1), rmax(2));
        pos[1] = Vector_t(0.0, 0.0, sposRef);
        pos[2] = Vector_t(rmin(0), rmin(1), rmin(2));

        for(int k = 0; k < 3; ++k) {
            externalB = Vector_t(0.0);
            externalE = Vector_t(0.0);
            FDext[2*k]   = externalB;
            FDext[2*k+1] = externalE * 1e-6;
        }

341
        itsDataSink_m->writeStatData(*itsBunch_m, FDext, rmax(2), sposRef, rmin(2));
gsell's avatar
gsell committed
342
        INFOMSG("* Wrote beam statistics." << endl);
343
    } else
gsell's avatar
gsell committed
344 345 346
        INFOMSG("* Invalid bunch! No statistics dumped." << endl);
}

347
#endif