ParallelSliceTracker.h 10.9 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>

kraus's avatar
kraus committed
7
#include "Algorithms/bet/EnvelopeBunch.h"
8

9
#include "Algorithms/Vektor.h"
gsell's avatar
gsell committed
10 11 12
#include "Algorithms/Tracker.h"
#include "Structure/DataSink.h"
#include "Utilities/Options.h"
13
#include "Utilities/Options.h"
gsell's avatar
gsell committed
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 56

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

    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
68
    virtual void visitDegrader(const Degrader &);
gsell's avatar
gsell committed
69 70 71 72 73
    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 &);
74
    virtual void visitMultipole(const Multipole &);
gsell's avatar
gsell committed
75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92
    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 &);

93

gsell's avatar
gsell committed
94 95 96 97 98 99
private:

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

100
    void updateRFElement(std::string elName, double maxPhi);
101
    void printRFPhases();
gsell's avatar
gsell committed
102

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
    DataSink *itsDataSink_m;
gsell's avatar
gsell committed
111 112

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

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

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

    void kickParticles();

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

135
    void prepareSections();
136
    void handleAutoPhasing();
137 138 139 140 141 142 143 144 145 146 147
    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
148 149 150 151 152 153 154 155 156 157 158 159
    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) {
160
    itsOpalBeamline_m->visit(wrap, *this, itsBunch_m);
gsell's avatar
gsell committed
161 162 163
}

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

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


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

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

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


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


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


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


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


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


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


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


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

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


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

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


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


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


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

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() {
264
    itsBunch_m->calcBeamParameters();
gsell's avatar
gsell committed
265 266 267 268

}

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

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;
278
    Vector_t FDext[6];  // = {BHead, EHead, BRef, ERef, BTail, ETail}.
gsell's avatar
gsell committed
279

280 281 282 283 284
    /*
     * 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
285
    Vector_t rmin, rmax;
286
    itsBunch_m->get_bounds(rmin, rmax);
gsell's avatar
gsell committed
287 288 289 290 291 292 293 294 295

    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);
296 297 298
        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
299 300 301 302 303 304

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

    if(step % Options::psDumpFreq == 0) {
305 306 307 308
        //itsDataSink->stashPhaseSpaceEnvelope(*itsBunch_m, FDext, rmax(2), sposRef, rmin(2));
        itsDataSink_m->writePhaseSpaceEnvelope(*itsBunch_m, FDext,
                                             rmax(2), sposRef, rmin(2));
        msg     << *itsBunch_m << endl;
gsell's avatar
gsell committed
309 310 311
    }

    if(step % Options::statDumpFreq == 0) {
312 313
        itsDataSink_m->writeStatData(*itsBunch_m, FDext,
                                   rmax(2), sposRef, rmin(2));
gsell's avatar
gsell committed
314 315 316 317 318
    }
}

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

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

337 338
        itsDataSink_m->writeStatData(*itsBunch_m, FDext, rmax(2), sposRef, rmin(2));
    } else
gsell's avatar
gsell committed
339 340 341
        INFOMSG("* Invalid bunch! No statistics dumped." << endl);
}

342
#endif