ThickTracker.h

Go to the documentation of this file.

#ifndef OPAL_ThickTracker_HH
#define OPAL_ThickTracker_HH

// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: ThickTracker
//
// ------------------------------------------------------------------------
//
// $Author: ganz_p $
//
// ------------------------------------------------------------------------

#include "Algorithms/Tracker.h"
#include "Structure/DataSink.h"

#include "Hamiltonian.h"

#include "MapAnalyser.h"

#include "Algorithms/IndexMap.h"
#include "AbsBeamline/BeamBeam.h"
#include "AbsBeamline/BeamStripping.h"
#include "AbsBeamline/CCollimator.h"
#include "AbsBeamline/Corrector.h"
#include "AbsBeamline/CyclotronValley.h"
#include "AbsBeamline/Diagnostic.h"
#include "AbsBeamline/Degrader.h"
#include "AbsBeamline/Drift.h"
#include "AbsBeamline/FlexibleCollimator.h"
#include "AbsBeamline/ElementBase.h"
#include "AbsBeamline/Lambertson.h"
#include "AbsBeamline/Marker.h"
#include "AbsBeamline/Monitor.h"
#include "AbsBeamline/Multipole.h"
#include "AbsBeamline/ParallelPlate.h"
#include "AbsBeamline/Probe.h"
#include "AbsBeamline/RFCavity.h"
#include "AbsBeamline/RFQuadrupole.h"
#include "AbsBeamline/RBend.h"
#include "AbsBeamline/RBend3D.h"
#include "AbsBeamline/SBend.h"
#include "AbsBeamline/Separator.h"
#include "AbsBeamline/Septum.h"
#include "AbsBeamline/Solenoid.h"
#include "AbsBeamline/TravelingWave.h"

#include "Elements/OpalBeamline.h"

#include "Structure/Beam.h"

//#include <array>
#include <cmath>

#include <tuple>

class BMultipoleField;

template <class T, unsigned Dim>
class PartBunchBase;



// Class ThickTracker
// ------------------------------------------------------------------------
// [p]
// Phase space coordinates numbering:
// [tab 3 b]
// [row]number [&]name          [&]unit  [/row]
// [row]0      [&]$x$           [&]metres [/row]
// [row]1      [&]$p_x/p_r$     [&]1      [/row]
// [row]2      [&]$y$           [&]metres [/row]
// [row]3      [&]$p_y/p_r$     [&]1      [/row]
// [row]4      [&]$v*delta_t$   [&]metres [/row]
// [row]5      [&]$delta_p/p_r$ [&]1      [/row]
// [/tab][p]
// Where $p_r$ is the constant reference momentum defining the reference
// frame velocity, $m$ is the rest mass of the particles, and $v$ is the
// instantaneous velocity of the particle.
// [p]
// Other units used:
// [tab 2 b]
// [row]quantity             [&]unit           [/row]
// [row]reference momentum   [&]electron-volts [/row]
// [row]velocity             [&]metres/second  [/row]
// [row]accelerating voltage [&]volts          [/row]
// [row]separator voltage    [&]volts          [/row]
// [row]frequencies          [&]hertz          [/row]
// [row]phase lags           [&]$2*pi$         [/row]
// [/tab][p]
// Approximations used:
// [ul]
// [li]
// [li]
// [li]
// [/ul]
//
// On going through an element, we use the following steps:
// To complete the map, we propagate the closed orbit and add that to the map.

class ThickTracker: public Tracker {

public:
    typedef Hamiltonian::series_t                       series_t;
    typedef MapAnalyser::fMatrix_t                      fMatrix_t;
    typedef MapAnalyser::cfMatrix_t                     cfMatrix_t;
    typedef FVps<double, 6>                             map_t;
    typedef FVector<double, 6>                          particle_t;
    typedef std::tuple<series_t, std::size_t, double>   tuple_t;
    typedef std::list<tuple_t>                          beamline_t;

public:

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

    //  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 ThickTracker(const Beamline &bl,
                          PartBunchBase<double, 3> *bunch,
                          Beam &beam,
                          DataSink &ds,
                          const PartData &data,
                          bool revBeam, bool revTrack,
                          const std::vector<unsigned long long> &maxSTEPS,
                          double zstart,
                          const std::vector<double> &zstop,
                          const std::vector<double> &dt,
                          const int& truncOrder);

    virtual ~ThickTracker();


    virtual void visitAlignWrapper(const AlignWrapper &);

    virtual void visitBeamBeam(const BeamBeam &);

    virtual void visitBeamStripping(const BeamStripping &);

    virtual void visitCCollimator(const CCollimator &);

    virtual void visitCorrector(const Corrector &);

    virtual void visitDegrader(const Degrader &);

    virtual void visitDiagnostic(const Diagnostic &);

    virtual void visitDrift(const Drift &);

    virtual void visitFlexibleCollimator(const FlexibleCollimator &);

    virtual void visitLambertson(const Lambertson &);

    virtual void visitMarker(const Marker &);

    virtual void visitMonitor(const Monitor &);

    virtual void visitMultipole(const Multipole &);
    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 &);



    //  overwrite the execute-methode from DefaultVisitor
    virtual void visitBeamline(const Beamline &);

    //  overwrite the execute-methode from DefaultVisitor
    virtual void execute();

    void prepareSections();

    /*
    void insertFringeField(SBend* pSBend, std::list<structMapTracking>& mBL, double& beta0,
             double& gamma0, double& P0, double& q, std::array<double,2>& entrFringe, std::string e);
     */
private:


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

    void throwElementError_m(std::string element) {
        throw LogicalError("ThickTracker::execute()",
                           "Element '" + element + "' not supported.");
    }



    void checkElementOrder_m();

    void fillGaps_m();

    void track_m();

    particle_t particleToVector_m(const Vector_t& R,
                                  const Vector_t& P) const;

    void vectorToParticle_m(const particle_t& particle,
                            Vector_t& R,
                            Vector_t& P) const;

    void advanceParticles_m(const map_t& map);

    void updateParticle_m(particle_t& particle,
                          const map_t& map);

    void dump_m();


    void update_m(const double& spos,
                  const std::size_t& step);


    void write_m(const map_t& map);

    void concatenateMaps_m(const map_t& x, map_t& y);


    void advanceDispersion_m(fMatrix_t tempMatrix,
                             FMatrix<double, 1, 4> initialVal,
                             double pos);

    void applyEntranceFringe(double edge, double curve,
                             const BMultipoleField &field, double scale);

    void applyExitFringe(double edge, double curve,
                         const BMultipoleField &field, double scale);


    Hamiltonian hamiltonian_m;
    MapAnalyser mapAnalyser_m;


    Vector_t RefPartR_m;
    Vector_t RefPartP_m;

    DataSink *itsDataSink_m;

    OpalBeamline itsOpalBeamline_m;

    double zstart_m;        
    double zstop_m;         
    double threshold_m;     
    beamline_t elements_m;  

    CoordinateSystemTrafo referenceToLabCSTrafo_m;


    int truncOrder_m; 

    IpplTimings::TimerRef mapCreation_m;    
    IpplTimings::TimerRef mapCombination_m; 
    IpplTimings::TimerRef mapTracking_m;    
};

inline void ThickTracker::visitAlignWrapper(const AlignWrapper &wrap) {
//     itsOpalBeamline_m.visit(wrap, *this, itsBunch_m);
    this->throwElementError_m("AlignWrapper");
}

inline void ThickTracker::visitBeamBeam(const BeamBeam &bb) {
//     itsOpalBeamline_m.visit(bb, *this, itsBunch_m);
    this->throwElementError_m("BeamBeam");
}

inline void ThickTracker::visitBeamStripping(const BeamStripping &bstp) {
    itsOpalBeamline_m.visit(bstp, *this, itsBunch_m);
}

inline void ThickTracker::visitCCollimator(const CCollimator &coll) {
//     itsOpalBeamline_m.visit(coll, *this, itsBunch_m);
    this->throwElementError_m("CCollimator");
}


inline void ThickTracker::visitCorrector(const Corrector &corr) {
//     itsOpalBeamline_m.visit(corr, *this, itsBunch_m);
    this->throwElementError_m("Corrector");
}


inline void ThickTracker::visitDegrader(const Degrader &deg) {
//     itsOpalBeamline_m.visit(deg, *this, itsBunch_m);
    this->throwElementError_m("Degrader");
}


inline void ThickTracker::visitDiagnostic(const Diagnostic &diag) {
//     itsOpalBeamline_m.visit(diag, *this, itsBunch_m);
    this->throwElementError_m("Diagnostic");
}


inline void ThickTracker::visitDrift(const Drift &drift) {
    itsOpalBeamline_m.visit(drift, *this, itsBunch_m);


    double gamma = itsReference.getGamma();
    std::size_t nSlices = drift.getNSlices();
    double length       = drift.getElementLength();

    elements_m.push_back(std::make_tuple(hamiltonian_m.drift(gamma),
                                         nSlices,
                                         length));
}


inline void ThickTracker::visitFlexibleCollimator(const FlexibleCollimator &coll) {
//     itsOpalBeamline_m.visit(coll, *this, itsBunch_m);
    this->throwElementError_m("FlexibleCollimator");
}


inline void ThickTracker::visitLambertson(const Lambertson &lamb) {
//     itsOpalBeamline_m.visit(lamb, *this, itsBunch_m);
    this->throwElementError_m("Lambertson");
}


inline void ThickTracker::visitMarker(const Marker &marker) {
//     itsOpalBeamline_m.visit(marker, *this, itsBunch_m);
//     this->throwElementError_m("Marker");
}


inline void ThickTracker::visitMonitor(const Monitor &mon) {
//     itsOpalBeamline_m.visit(mon, *this, itsBunch_m);
    this->throwElementError_m("Monitor");
}


inline void ThickTracker::visitMultipole(const Multipole &mult) {
    itsOpalBeamline_m.visit(mult, *this, itsBunch_m);

    std::size_t nSlices = mult.getNSlices();
    double length       = mult.getElementLength();
    double gamma        = itsReference.getGamma();
    double p0           = itsReference.getP();
    double q            = itsReference.getQ(); // particle change [e]
    double gradB        = mult.getField().getNormalComponent(2) * ( Physics::c/ p0 ); // [T / m]
    //FIXME remove the next line
    gradB = std::round(gradB*1e6)*1e-6;

    double k1           = gradB * q *Physics::c / p0; // [1 / m^2]

    elements_m.push_back(std::make_tuple(hamiltonian_m.quadrupole(gamma, q, k1),
                                         nSlices,
                                         length));
}

inline void ThickTracker::visitProbe(const Probe &prob) {
//     itsOpalBeamline_m.visit(prob, *this, itsBunch_m);
    this->throwElementError_m("Probe");
}


inline void ThickTracker::visitRBend(const RBend &bend) {
//     itsOpalBeamline_m.visit(bend, *this, itsBunch_m);
    this->throwElementError_m("RBend");
}


inline void ThickTracker::visitRFCavity(const RFCavity &as) {
//     itsOpalBeamline_m.visit(as, *this, itsBunch_m);
    this->throwElementError_m("RFCavity");
}

inline void ThickTracker::visitTravelingWave(const TravelingWave &as) {
//     itsOpalBeamline_m.visit(as, *this, itsBunch_m);
    this->throwElementError_m("TravelingWave");
}


inline void ThickTracker::visitRFQuadrupole(const RFQuadrupole &rfq) {
//     itsOpalBeamline_m.visit(rfq, *this, itsBunch_m);
    this->throwElementError_m("RFQuadrupole");
}

inline void ThickTracker::visitSBend(const SBend &bend) {
    itsOpalBeamline_m.visit(bend, *this, itsBunch_m);

    double q     = itsReference.getQ(); // particle change [e]
    double ekin = bend.getDesignEnergy();

    double m     = itsReference.getM(); // eV / c^2
    double gamma = ekin / m + 1.0;
    double beta  = std::sqrt(1.0 - 1.0 / ( gamma * gamma ) );
    double p0    = gamma * beta * m; // eV / c

    double B     = bend.getB() * Physics::c / p0; // T = V * s / m^2
    double r     = std::abs(p0   / ( q * B * Physics::c )); // m

    double k0    = B * q * Physics::c / p0; // V * s * e * m / (m^2 * s * c )

    // [1/m]
    double h = 1.0 / r;

    double L = bend.getElementLength();

    if ( k0 < 0.0 )
        h *= -1.0;

    std::size_t nSlices = bend.getNSlices();

    double arclength    = 2.0 * r * std::asin( L / ( 2.0 * r ) ); //bend.getEffectiveLength();

    // Fringe fields currently not working
    //FIXME e1 not initialised
    //insert Entrance Fringefield
    double e1 = bend.getEntranceAngle();
    if (std::abs(e1) > 1e-6){
        elements_m.push_back(std::make_tuple(hamiltonian_m.fringeField(e1, h),
                                             1, 0.0));
    }

    //insert Dipole "body"
    elements_m.push_back(std::make_tuple(hamiltonian_m.sbend(gamma, h, k0),
                                         nSlices,
                                         arclength));

    //FIXME e2 not initialised
    //insert Exit Fringe field
    double e2 = bend.getExitAngle();
    if (std::abs(e2) > 1e-6){
        elements_m.push_back(std::make_tuple(hamiltonian_m.fringeField(e2, h),
                                             1, 0.0));
    }

}


inline void ThickTracker::visitSeparator(const Separator &sep) {
//     itsOpalBeamline_m.visit(sep, *this, itsBunch_m);
    this->throwElementError_m("Separator");
}


inline void ThickTracker::visitSeptum(const Septum &sept) {
//     itsOpalBeamline_m.visit(sept, *this, itsBunch_m);
    this->throwElementError_m("Septum");
}


inline void ThickTracker::visitSolenoid(const Solenoid &solenoid) {
//     itsOpalBeamline_m.visit(solenoid, *this, itsBunch_m);
    this->throwElementError_m("Solenoid");
}


inline void ThickTracker::visitParallelPlate(const ParallelPlate &pplate) {
//     itsOpalBeamline_m.visit(pplate, *this, itsBunch_m);
    this->throwElementError_m("ParallelPlate");
}

inline void ThickTracker::visitCyclotronValley(const CyclotronValley &cv) {
//     itsOpalBeamline_m.visit(cv, *this, itsBunch_m);
    this->throwElementError_m("CyclotronValley");
}

#endif // OPAL_ThickTracker_HH