4
#include "CLHEP/Matrix/Matrix.h"
5
#include "CLHEP/Matrix/Vector.h"
6
#include "CLHEP/Units/SystemOfUnits.h"
8
#include "src/common/Interface/Differentiator.hh"
9
#include "CovarianceMatrix.hh"
10
#include "PhaseSpaceVector.hh"
11
#include "OpticsModel.hh"
15
//TransferMap is a polynomial mapping of a phase space vector from a plane at Z1 to another plane at Z2,
16
//such that PhaseSpaceVector with coordinates U transform like
17
// U(Z2)_p = M_ip U(Z1)_i + M_ijp U(Z1)_i U(Z1)_j + M_ijkp U(Z1)_i U(Z1)_j U(Z1)_k + ...
18
//extending to arbitrary order. A reference trajectory is assumed, that is the transformation is applied about
19
//some phase space vector that is taken to be a zero point.
21
//Currently first order transformations use the CLHEP::HepMatrix object, second order transformations use the
22
//Tensor3 object and higher order transformations use the Tensor object
24
//At some point I would like to move to the more general PolynomialVector object. For now I have implemented an
25
//interface to this object
32
//Create the transfer map that transports from referenceIn to referenceOut using firstOrderMap
33
TransferMap(HepMatrix firstOrderMap, PhaseSpaceVector referenceIn, PhaseSpaceVector referenceOut,
35
TransferMap(HepMatrix transverseMap, HepMatrix longitudinalMap, PhaseSpaceVector referenceIn,
36
PhaseSpaceVector referenceOut, OpticsModel* Optics);
37
//Create a transfer map that does nothing i.e. referenceIn=referenceOut, M=1
38
TransferMap(PhaseSpaceVector referenceIn, OpticsModel* Optics);
42
void SetFirstOrderMap(const HepMatrix transverseMap, const HepMatrix longitudinalMap);
43
void SetFirstOrderMap(const HepMatrix firstOrderMap);
44
void SetSecondOrderMap(Tensor3 secondOrderMap) {_secondOrderMap = secondOrderMap;}
45
void SetReferenceIn(const PhaseSpaceVector referenceIn) {_referenceTrajectoryIn = referenceIn;}
46
void SetReferenceOut(const PhaseSpaceVector referenceOut) {_referenceTrajectoryOut = referenceOut;}
47
void SetReferenceOutZ(double z) {_referenceTrajectoryOut.setZ(z);}
48
void IncrementReferenceOutZ(double dz) {SetReferenceOutZ(_referenceTrajectoryOut.z()+dz);}
49
void SetOpticsModel(OpticsModel* optics) {_optics = optics;}
50
//Accessors or whatever
51
PhaseSpaceVector GetReferenceIn() const {return _referenceTrajectoryIn;}
52
PhaseSpaceVector GetReferenceOut() const {return _referenceTrajectoryOut;}
53
//Accessor for first order matrix
54
HepMatrix GetFirstOrderMap() const;
55
HepMatrix GetTransverseMap() const {return _firstOrderMap.sub(3,6,3,6);}
56
HepMatrix GetLongitudinalMap() const {return _firstOrderMap.sub(1,2,1,2);}
57
//Accessor for second and higher order maps
58
Tensor3 GetSecondOrderMap() const {return _secondOrderMap;}
59
Tensor GetMap(int order) const;
60
//Methods to transport objects
61
//Return value is M^T V M
62
CovarianceMatrix operator*(const CovarianceMatrix& aCovMatrix) const;
64
PhaseSpaceVector operator*(const PhaseSpaceVector& aPhaseSpaceVector) const;
65
//Return value is M1*M2 where *this is M1 and aTransferMAp is M2
66
TransferMap operator*(const TransferMap& aTransferMap) const;
67
TransferMap operator*=(const TransferMap& aTransferMap);
68
//Return transfer map rotated by an angle theta
69
TransferMap Rotate(double angle) const;
70
//Return the inverse of the transfer map
71
TransferMap Inv(int & failflag) const;
72
TransferMap Inv() const {int i; return Inv(i);}
73
//determinant of the 6d transfer map
74
double Det() const {return GetFirstOrderMap().determinant();}
75
//decompose into a "matched" covariance matrix
76
CovarianceMatrix Decompose(double em_t, double em_x, double em_y) const;
77
//rotate by angle then decompose into a "matched" covariance matrix... then rotate cov matrix back
78
CovarianceMatrix Decompose(double angle, double em_t, double em_x, double em_y) const;
79
//return a cov matrix with coefficients for sxx_out/sxx_in, spp_out/spp_in and sxp_out/sxp_in
80
std::vector<CovarianceMatrix> Decompose(double sxx, double sxp, double spp, double em_t, double em_x,
82
CovarianceMatrix Decompose(double sxx, double sxp, double spp,
83
double em_t, double em_x, double em_y, double angle) const;
84
//Get the phase advance of the transfer matrix in the domain (0, pi), given by acos((m[0][0]+m[1][1])/2.);
85
//units are radians; if LarmorAngle is specified, unrotates matrix first
86
double PhaseAdvance(int axis) const;
87
double PhaseAdvance(int axis, double LarmorAngle) const;
89
PolynomialVector* GetPolynomialVector() {return _polynomial;}
90
void SetPolynomialVector(PolynomialVector* poly) {_polynomial = poly;}
91
bool IsCanonical() {return _canonicalMap;}
92
bool IsCanonical(bool can) {_canonicalMap = can; return _canonicalMap;}
95
static void SetOrder(int order) {_order = order;}
96
static int GetOrder() {return _order;}
97
friend std::ostream& operator<<(std::ostream& out, TransferMap tm);
100
//Transfer maps of various orders
101
HepMatrix _firstOrderMap;
102
Tensor3 _secondOrderMap;
103
std::vector<Tensor*> _higherOrderMaps;
104
//The reference trajectory is a property of the transfer map
105
PhaseSpaceVector _referenceTrajectoryIn;
106
PhaseSpaceVector _referenceTrajectoryOut;
107
//Need this for transporting phase space vectors
108
OpticsModel* _optics;
109
PolynomialVector* _polynomial; //set to NULL to disable
115
std::ostream& operator<<(std::ostream& out, TransferMap tm);
added
removed
tests/integration/optics/src/TransferMap.hh
