Skip to content
GitLab
Projects
Groups
Snippets
Help
Loading...
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Sign in
Toggle navigation
Open sidebar
OPAL
src
Commits
d9203ef7
Commit
d9203ef7
authored
Jul 14, 2020
by
frey_m
Browse files
generate matched distribution with symplectic method
parent
b780b391
Changes
8
Expand all
Hide whitespace changes
Inline
Side-by-side
Showing
8 changed files
with
714 additions
and
545 deletions
+714
-545
src/Distribution/CMakeLists.txt
src/Distribution/CMakeLists.txt
+2
-2
src/Distribution/Distribution.cpp
src/Distribution/Distribution.cpp
+142
-52
src/Distribution/Distribution.h
src/Distribution/Distribution.h
+4
-0
src/Distribution/RDM.h
src/Distribution/RDM.h
+0
-99
src/Distribution/RealDiracMatrix.cpp
src/Distribution/RealDiracMatrix.cpp
+343
-0
src/Distribution/RealDiracMatrix.h
src/Distribution/RealDiracMatrix.h
+123
-0
src/Distribution/SigmaGenerator.cpp
src/Distribution/SigmaGenerator.cpp
+63
-318
src/Distribution/SigmaGenerator.h
src/Distribution/SigmaGenerator.h
+37
-74
No files found.
src/Distribution/CMakeLists.txt
View file @
d9203ef7
set
(
_SRCS
Distribution.cpp
LaserProfile.cpp
R
DM
.cpp
R
ealDiracMatrix
.cpp
SigmaGenerator.cpp
)
...
...
@@ -17,7 +17,7 @@ set (HDRS
LaserProfile.h
MapGenerator.h
matrix_vector_operation.h
R
DM
.h
R
ealDiracMatrix
.h
SigmaGenerator.h
)
...
...
src/Distribution/Distribution.cpp
View file @
d9203ef7
...
...
@@ -7,7 +7,6 @@
// OPAL is licensed under GNU GPL version 3.
#include "Distribution/Distribution.h"
#include "Distribution/SigmaGenerator.h"
#include "Distribution/ClosedOrbitFinder.h"
#include "AbsBeamline/SpecificElementVisitor.h"
...
...
@@ -1302,58 +1301,18 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
for
(
unsigned
int
i
=
0
;
i
<
siggen
->
getSigma
().
size1
();
++
i
)
{
*
gmsg
<<
std
::
setprecision
(
4
)
<<
std
::
setw
(
11
)
<<
siggen
->
getSigma
()(
i
,
0
);
for
(
unsigned
int
j
=
1
;
j
<
siggen
->
getSigma
().
size2
();
++
j
)
{
if
(
siggen
->
getSigma
()(
i
,
j
)
<
10e-1
2
)
{
if
(
std
::
abs
(
siggen
->
getSigma
()(
i
,
j
)
)
<
1
.
0e-1
5
)
{
*
gmsg
<<
" & "
<<
std
::
setprecision
(
4
)
<<
std
::
setw
(
11
)
<<
0.0
;
}
else
{
*
gmsg
<<
" & "
<<
std
::
setprecision
(
4
)
<<
std
::
setw
(
11
)
<<
siggen
->
getSigma
()(
i
,
j
);
*
gmsg
<<
" & "
<<
std
::
setprecision
(
4
)
<<
std
::
setw
(
11
)
<<
siggen
->
getSigma
()(
i
,
j
);
}
}
*
gmsg
<<
"
\\\\
"
<<
endl
;
}
/*
Now setup the distribution generator
Units of the Sigma Matrix:
spatial: m
moment: rad
*/
double
gamma
=
E_m
/
massIneV
+
1.0
;
double
beta
=
std
::
sqrt
(
1.0
-
1.0
/
(
gamma
*
gamma
));
auto
sigma
=
siggen
->
getSigma
();
for
(
unsigned
int
i
=
0
;
i
<
3
;
++
i
)
{
if
(
sigma
(
2
*
i
,
2
*
i
)
<
0
||
sigma
(
2
*
i
+
1
,
2
*
i
+
1
)
<
0
)
throw
OpalException
(
"Distribution::CreateMatchedGaussDistribution()"
,
"Negative value on the diagonal of the sigma matrix."
);
}
sigmaR_m
[
0
]
=
std
::
sqrt
(
sigma
(
0
,
0
));
sigmaP_m
[
0
]
=
std
::
sqrt
(
sigma
(
1
,
1
))
*
beta
*
gamma
;
sigmaR_m
[
2
]
=
std
::
sqrt
(
sigma
(
2
,
2
));
sigmaP_m
[
2
]
=
std
::
sqrt
(
sigma
(
3
,
3
))
*
beta
*
gamma
;
sigmaR_m
[
1
]
=
std
::
sqrt
(
sigma
(
4
,
4
));
//p_l^2 = [(delta+1)*beta*gamma]^2 - px^2 - pz^2
double
pl2
=
(
std
::
sqrt
(
sigma
(
5
,
5
))
+
1
)
*
(
std
::
sqrt
(
sigma
(
5
,
5
))
+
1
)
*
beta
*
gamma
*
beta
*
gamma
-
sigmaP_m
[
0
]
*
sigmaP_m
[
0
]
-
sigmaP_m
[
2
]
*
sigmaP_m
[
2
];
double
pl
=
std
::
sqrt
(
pl2
);
sigmaP_m
[
1
]
=
gamma
*
(
pl
-
beta
*
gamma
);
correlationMatrix_m
(
1
,
0
)
=
sigma
(
0
,
1
)
/
(
std
::
sqrt
(
sigma
(
0
,
0
)
*
sigma
(
1
,
1
)));
correlationMatrix_m
(
3
,
2
)
=
sigma
(
2
,
3
)
/
(
std
::
sqrt
(
sigma
(
2
,
2
)
*
sigma
(
3
,
3
)));
correlationMatrix_m
(
5
,
4
)
=
sigma
(
4
,
5
)
/
(
std
::
sqrt
(
sigma
(
4
,
4
)
*
sigma
(
5
,
5
)));
correlationMatrix_m
(
4
,
0
)
=
sigma
(
0
,
4
)
/
(
std
::
sqrt
(
sigma
(
0
,
0
)
*
sigma
(
4
,
4
)));
correlationMatrix_m
(
4
,
1
)
=
sigma
(
1
,
4
)
/
(
std
::
sqrt
(
sigma
(
1
,
1
)
*
sigma
(
4
,
4
)));
correlationMatrix_m
(
5
,
0
)
=
sigma
(
0
,
5
)
/
(
std
::
sqrt
(
sigma
(
0
,
0
)
*
sigma
(
5
,
5
)));
correlationMatrix_m
(
5
,
1
)
=
sigma
(
1
,
5
)
/
(
std
::
sqrt
(
sigma
(
1
,
1
)
*
sigma
(
5
,
5
)));
createDistributionGauss
(
numberOfParticles
,
massIneV
);
generateMatchedGauss
(
siggen
->
getSigma
(),
numberOfParticles
,
massIneV
);
// update injection radius and radial momentum
CyclotronElement
->
setRinit
(
siggen
->
getInjectionRadius
()
*
1.0e3
);
...
...
@@ -2404,6 +2363,137 @@ void Distribution::generateGaussZ(size_t numberOfParticles) {
gsl_matrix_free
(
corMat
);
}
void
Distribution
::
generateMatchedGauss
(
const
SigmaGenerator
::
matrix_t
&
sigma
,
size_t
numberOfParticles
,
double
massIneV
)
{
/* This particle distribution generation is based on a
* symplectic method described in
* https://arxiv.org/abs/1205.3601
*/
/* Units of the Sigma Matrix:
* spatial: m
* moment: rad
*
* Attention: The vertical and longitudinal direction must be interchanged!
*/
for
(
unsigned
int
i
=
0
;
i
<
3
;
++
i
)
{
if
(
sigma
(
2
*
i
,
2
*
i
)
<
0
||
sigma
(
2
*
i
+
1
,
2
*
i
+
1
)
<
0
)
throw
OpalException
(
"Distribution::generateMatchedGauss()"
,
"Negative value on the diagonal of the sigma matrix."
);
}
double
gamma
=
E_m
/
massIneV
+
1.0
;
// beta * gamma
double
bgam
=
std
::
sqrt
(
gamma
*
gamma
-
1.0
);
/*
* only used for printing
*/
// horitzonal
sigmaR_m
[
0
]
=
std
::
sqrt
(
sigma
(
0
,
0
));
sigmaP_m
[
0
]
=
std
::
sqrt
(
sigma
(
1
,
1
))
*
bgam
;
// longitudinal
sigmaR_m
[
1
]
=
std
::
sqrt
(
sigma
(
4
,
4
));
sigmaP_m
[
1
]
=
std
::
sqrt
(
sigma
(
5
,
5
))
*
bgam
;
// vertical
sigmaR_m
[
2
]
=
std
::
sqrt
(
sigma
(
2
,
2
));
sigmaP_m
[
2
]
=
std
::
sqrt
(
sigma
(
3
,
3
))
*
bgam
;
correlationMatrix_m
(
1
,
0
)
=
sigma
(
0
,
1
)
/
(
std
::
sqrt
(
sigma
(
0
,
0
)
*
sigma
(
1
,
1
)));
correlationMatrix_m
(
3
,
2
)
=
sigma
(
2
,
3
)
/
(
std
::
sqrt
(
sigma
(
2
,
2
)
*
sigma
(
3
,
3
)));
correlationMatrix_m
(
5
,
4
)
=
sigma
(
4
,
5
)
/
(
std
::
sqrt
(
sigma
(
4
,
4
)
*
sigma
(
5
,
5
)));
correlationMatrix_m
(
4
,
0
)
=
sigma
(
0
,
4
)
/
(
std
::
sqrt
(
sigma
(
0
,
0
)
*
sigma
(
4
,
4
)));
correlationMatrix_m
(
4
,
1
)
=
sigma
(
1
,
4
)
/
(
std
::
sqrt
(
sigma
(
1
,
1
)
*
sigma
(
4
,
4
)));
correlationMatrix_m
(
5
,
0
)
=
sigma
(
0
,
5
)
/
(
std
::
sqrt
(
sigma
(
0
,
0
)
*
sigma
(
5
,
5
)));
correlationMatrix_m
(
5
,
1
)
=
sigma
(
1
,
5
)
/
(
std
::
sqrt
(
sigma
(
1
,
1
)
*
sigma
(
5
,
5
)));
inputMoUnits_m
=
InputMomentumUnitsT
::
NONE
;
/*
* decouple horitzonal and longitudinal direction
*/
// extract horitzonal and longitudinal directions
RealDiracMatrix
::
matrix_t
A
(
4
,
4
);
A
(
0
,
0
)
=
sigma
(
0
,
0
);
A
(
1
,
1
)
=
sigma
(
1
,
1
);
A
(
2
,
2
)
=
sigma
(
4
,
4
);
A
(
3
,
3
)
=
sigma
(
5
,
5
);
A
(
0
,
1
)
=
sigma
(
0
,
1
);
A
(
0
,
2
)
=
sigma
(
0
,
4
);
A
(
0
,
3
)
=
sigma
(
0
,
5
);
A
(
1
,
0
)
=
sigma
(
1
,
0
);
A
(
2
,
0
)
=
sigma
(
4
,
0
);
A
(
3
,
0
)
=
sigma
(
5
,
0
);
A
(
1
,
2
)
=
sigma
(
1
,
4
);
A
(
2
,
1
)
=
sigma
(
4
,
1
);
A
(
1
,
3
)
=
sigma
(
1
,
5
);
A
(
3
,
1
)
=
sigma
(
5
,
1
);
A
(
2
,
3
)
=
sigma
(
4
,
5
);
A
(
3
,
2
)
=
sigma
(
5
,
4
);
RealDiracMatrix
rdm
;
RealDiracMatrix
::
sparse_matrix_t
R1
=
rdm
.
diagonalize
(
A
);
RealDiracMatrix
::
vector_t
variances
(
8
);
for
(
int
i
=
0
;
i
<
4
;
++
i
)
{
variances
(
i
)
=
std
::
sqrt
(
A
(
i
,
i
));
}
/*
* decouple vertical direction
*/
A
*=
0.0
;
A
(
0
,
0
)
=
1
;
A
(
1
,
1
)
=
1
;
A
(
2
,
2
)
=
sigma
(
2
,
2
);
A
(
3
,
3
)
=
sigma
(
3
,
3
);
A
(
2
,
3
)
=
sigma
(
2
,
3
);
A
(
3
,
2
)
=
sigma
(
3
,
2
);
RealDiracMatrix
::
sparse_matrix_t
R2
=
rdm
.
diagonalize
(
A
);
for
(
int
i
=
0
;
i
<
4
;
++
i
)
{
variances
(
4
+
i
)
=
std
::
sqrt
(
A
(
i
,
i
));
}
xDist_m
.
resize
(
numberOfParticles
);
pxDist_m
.
resize
(
numberOfParticles
);
yDist_m
.
resize
(
numberOfParticles
);
pyDist_m
.
resize
(
numberOfParticles
);
tOrZDist_m
.
resize
(
numberOfParticles
);
pzDist_m
.
resize
(
numberOfParticles
);
RealDiracMatrix
::
vector_t
p1
(
4
),
p2
(
4
);
for
(
size_t
i
=
0
;
i
<
numberOfParticles
;
i
++
)
{
for
(
int
j
=
0
;
j
<
4
;
j
++
)
{
p1
(
j
)
=
gsl_ran_gaussian
(
randGen_m
,
1.0
)
*
variances
(
j
);
p2
(
j
)
=
gsl_ran_gaussian
(
randGen_m
,
1.0
)
*
variances
(
4
+
j
);
}
p1
=
boost
::
numeric
::
ublas
::
prod
(
R1
,
p1
);
p2
=
boost
::
numeric
::
ublas
::
prod
(
R2
,
p2
);
xDist_m
[
i
]
=
p1
(
0
);
pxDist_m
[
i
]
=
p1
(
1
)
*
bgam
;
yDist_m
[
i
]
=
p1
(
2
);
pyDist_m
[
i
]
=
p1
(
3
)
*
bgam
;
tOrZDist_m
[
i
]
=
p2
(
2
);
pzDist_m
[
i
]
=
p2
(
3
)
*
bgam
;
}
}
void
Distribution
::
generateLongFlattopT
(
size_t
numberOfParticles
)
{
double
flattopTime
=
tPulseLengthFWHM_m
...
...
@@ -2963,7 +3053,7 @@ void Distribution::printDistMatchedGauss(Inform &os) const {
os
<<
"* SIGMAPX = "
<<
sigmaP_m
[
0
]
<<
" [Beta Gamma]"
<<
endl
;
os
<<
"* SIGMAPY = "
<<
sigmaP_m
[
1
]
<<
" [Beta Gamma]"
<<
endl
;
os
<<
"* SIGMAPZ = "
<<
sigmaP_m
[
2
]
<<
" [Beta Gamma]"
<<
endl
;
os
<<
"* AVRGPZ = "
<<
avrgpz_m
<<
" [Beta Gamma]"
<<
endl
;
//
os << "* AVRGPZ = " << avrgpz_m << " [Beta Gamma]" << endl;
os
<<
"* CORRX = "
<<
correlationMatrix_m
(
1
,
0
)
<<
endl
;
os
<<
"* CORRY = "
<<
correlationMatrix_m
(
3
,
2
)
<<
endl
;
...
...
@@ -2972,12 +3062,12 @@ void Distribution::printDistMatchedGauss(Inform &os) const {
os
<<
"* R62 = "
<<
correlationMatrix_m
(
5
,
1
)
<<
endl
;
os
<<
"* R51 = "
<<
correlationMatrix_m
(
4
,
0
)
<<
endl
;
os
<<
"* R52 = "
<<
correlationMatrix_m
(
4
,
1
)
<<
endl
;
os
<<
"* CUTOFFX = "
<<
cutoffR_m
[
0
]
<<
" [units of SIGMAX]"
<<
endl
;
os
<<
"* CUTOFFY = "
<<
cutoffR_m
[
1
]
<<
" [units of SIGMAY]"
<<
endl
;
os
<<
"* CUTOFFLONG = "
<<
cutoffR_m
[
2
]
<<
" [units of SIGMAZ]"
<<
endl
;
os
<<
"* CUTOFFPX = "
<<
cutoffP_m
[
0
]
<<
" [units of SIGMAPX]"
<<
endl
;
os
<<
"* CUTOFFPY = "
<<
cutoffP_m
[
1
]
<<
" [units of SIGMAPY]"
<<
endl
;
os
<<
"* CUTOFFPZ = "
<<
cutoffP_m
[
2
]
<<
" [units of SIGMAPY]"
<<
endl
;
//
os << "* CUTOFFX = " << cutoffR_m[0] << " [units of SIGMAX]" << endl;
//
os << "* CUTOFFY = " << cutoffR_m[1] << " [units of SIGMAY]" << endl;
//
os << "* CUTOFFLONG = " << cutoffR_m[2] << " [units of SIGMAZ]" << endl;
//
os << "* CUTOFFPX = " << cutoffP_m[0] << " [units of SIGMAPX]" << endl;
//
os << "* CUTOFFPY = " << cutoffP_m[1] << " [units of SIGMAPY]" << endl;
//
os << "* CUTOFFPZ = " << cutoffP_m[2] << " [units of SIGMAPY]" << endl;
}
void
Distribution
::
printDistGauss
(
Inform
&
os
)
const
{
...
...
@@ -3207,7 +3297,7 @@ void Distribution::setAttributes() {
// itsAttr[Attrib::Distribution::E2]
// = Attributes::makeReal("E2", "If E2<Eb, we compute the tunes from the beams energy Eb to E2 with dE=0.25 MeV ", 0.0);
itsAttr
[
Attrib
::
Distribution
::
RESIDUUM
]
=
Attributes
::
makeReal
(
"RESIDUUM"
,
"Residuum for the closed orbit finder and sigma matrix generator "
,
1e-
12
);
=
Attributes
::
makeReal
(
"RESIDUUM"
,
"Residuum for the closed orbit finder and sigma matrix generator "
,
1e-
8
);
itsAttr
[
Attrib
::
Distribution
::
MAXSTEPSCO
]
=
Attributes
::
makeReal
(
"MAXSTEPSCO"
,
"Maximum steps used to find closed orbit "
,
100
);
itsAttr
[
Attrib
::
Distribution
::
MAXSTEPSSI
]
...
...
src/Distribution/Distribution.h
View file @
d9203ef7
...
...
@@ -19,6 +19,8 @@
#include "AppTypes/SymTenzor.h"
#include "Attributes/Attributes.h"
#include "Distribution/SigmaGenerator.h"
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_qrng.h>
#include <gsl/gsl_rng.h>
...
...
@@ -336,6 +338,8 @@ private:
void
generateFlattopT
(
size_t
numberOfParticles
);
void
generateFlattopZ
(
size_t
numberOfParticles
);
void
generateGaussZ
(
size_t
numberOfParticles
);
void
generateMatchedGauss
(
const
SigmaGenerator
::
matrix_t
&
,
size_t
numberOfParticles
,
double
massIneV
);
void
generateLongFlattopT
(
size_t
numberOfParticles
);
void
generateTransverseGauss
(
size_t
numberOfParticles
);
void
initializeBeam
(
PartBunchBase
<
double
,
3
>
*
beam
);
...
...
src/Distribution/RDM.h
deleted
100644 → 0
View file @
b780b391
//
// Class: RDM
// They're ordered after the paper of Dr. C. Baumgarten: "Use of real Dirac matrices in two-dimensional coupled linear optics".
// The diagonalizing method is based on the paper "Geometrical method of decoupling" (2012) of Dr. C. Baumgarten.
//
// Copyright (c) 2014, Christian Baumgarten, Paul Scherrer Institut, Villigen PSI, Switzerland
// Matthias Frey, ETH Zürich
// All rights reserved
//
// Implemented as part of the Semester thesis by Matthias Frey
// "Matched Distributions in Cyclotrons"
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
#ifndef RDM_H
#define RDM_H
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.hpp>
/// @brief Real Dirac Matrix class
class
RDM
{
public:
/// Sparse matrix type definition
typedef
boost
::
numeric
::
ublas
::
compressed_matrix
<
double
,
boost
::
numeric
::
ublas
::
row_major
>
sparse_matrix_t
;
/// Dense matrix type definition
typedef
boost
::
numeric
::
ublas
::
matrix
<
double
>
matrix_t
;
/// Dense vector type definition
typedef
boost
::
numeric
::
ublas
::
vector
<
double
>
vector_t
;
/// Default constructor (sets only NumOfRDMs and DimOfRDMs)
RDM
();
/// Returns the i-th Real Dirac matrix
/*!
* @param i specifying the matrix (has to be in the range from 0 to 15)
*/
sparse_matrix_t
getRDM
(
short
);
/// Decomposes a real-valued 4x4 matrix into a linear combination and returns a vector containing the coefficients
/*!
* @param M an arbitrary real-valued 4x4 matrix
*/
vector_t
decompose
(
const
matrix_t
&
);
/// Takes a vector of coefficients, evaluates the linear combination of RDMs with these coefficients and returns a 4x4 matrix
/*!
* @param coeffs is a vector of coefficients (at most length NumOfRDMs)
*/
matrix_t
combine
(
const
vector_t
&
);
/// Brings a 4x4 symplex matrix into Hamilton form and computes the transformation matrix and its inverse
/*!
* @param Ms is a 4x4 symplex matrix
* @param R is the 4x4 transformation matrix (gets computed)
* @param invR is the 4x4 inverse transformation matrix (gets computed)
*/
void
diagonalize
(
matrix_t
&
,
sparse_matrix_t
&
,
sparse_matrix_t
&
);
/// Returns the symplex part of a 4x4 real-valued matrix
/*!
* @param M 4x4 real-valued matrix
*/
matrix_t
symplex
(
const
matrix_t
&
);
/// Returns the cosymplex part of a 4x4 real-valued matrix
/*!
* @param M 4x4 real-valued matrix
*/
matrix_t
cosymplex
(
const
matrix_t
&
);
/// The number of real Dirac matrices
short
NumOfRDMs
;
/// The matrix dimension (4x4)
short
DimOfRDMs
;
private:
/// Applies a rotation to the matrix M by a given angle
/*!
* @param M is the matrix to be transformed
* @param i is the i-th RDM used for transformation
* @param phi is the angle of rotation
* @param Rtot is a reference to the current transformation matrix
* @param invRtot is a reference to the inverse of the current transformation matrix
*/
void
transform
(
matrix_t
&
,
short
,
double
,
sparse_matrix_t
&
,
sparse_matrix_t
&
);
};
#endif
\ No newline at end of file
src/Distribution/R
DM
.cpp
→
src/Distribution/R
ealDiracMatrix
.cpp
View file @
d9203ef7
//
// Class: RDM
// Class: RealDiracMatrix
// Real Dirac matrix class.
// They're ordered after the paper of Dr. C. Baumgarten: "Use of real Dirac matrices in two-dimensional coupled linear optics".
// The diagonalizing method is based on the paper "Geometrical method of decoupling" (2012) of Dr. C. Baumgarten.
//
// Copyright (c) 2014, Christian Baumgarten, Paul Scherrer Institut, Villigen PSI, Switzerland
// Matthias Frey, ETH Zürich
// Copyright (c) 2014,
2020
Christian Baumgarten, Paul Scherrer Institut, Villigen PSI, Switzerland
//
Matthias Frey, ETH Zürich
// All rights reserved
//
// Implemented as part of the Semester thesis by Matthias Frey
...
...
@@ -20,36 +21,36 @@
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
#include "R
DM
.h"
#include "R
ealDiracMatrix
.h"
#include "Utilities/OpalException.h"
#include <cmath>
#include "matrix_vector_operation.h"
R
DM
::
RDM
()
:
NumOfRDMs
(
16
),
DimOfRDMs
(
4
)
{};
R
ealDiracMatrix
::
RealDiracMatrix
()
:
NumOfRDMs
(
16
),
DimOfRDMs
(
4
)
{};
typename
R
DM
::
sparse_matrix_t
R
DM
::
getRDM
(
short
i
)
{
typename
R
ealDiracMatrix
::
sparse_matrix_t
R
ealDiracMatrix
::
getRDM
(
short
i
)
{
sparse_matrix_t
rdm
(
4
,
4
,
4
);
// #nrows, #ncols, #non-zeros
switch
(
i
)
{
case
0
:
rdm
(
0
,
1
)
=
rdm
(
2
,
3
)
=
1
;
rdm
(
1
,
0
)
=
rdm
(
3
,
2
)
=
-
1
;
break
;
case
1
:
rdm
(
0
,
1
)
=
rdm
(
1
,
0
)
=
-
1
;
rdm
(
2
,
3
)
=
rdm
(
3
,
2
)
=
1
;
break
;
case
2
:
rdm
(
0
,
3
)
=
rdm
(
1
,
2
)
=
rdm
(
2
,
1
)
=
rdm
(
3
,
0
)
=
1
;
break
;
case
3
:
rdm
(
0
,
0
)
=
rdm
(
2
,
2
)
=
-
1
;
rdm
(
1
,
1
)
=
rdm
(
3
,
3
)
=
1
;
break
;
case
4
:
rdm
(
0
,
0
)
=
rdm
(
3
,
3
)
=
-
1
;
rdm
(
1
,
1
)
=
rdm
(
2
,
2
)
=
1
;
break
;
case
5
:
rdm
(
0
,
2
)
=
rdm
(
2
,
0
)
=
1
;
rdm
(
1
,
3
)
=
rdm
(
3
,
1
)
=
-
1
;
break
;
case
6
:
rdm
(
0
,
1
)
=
rdm
(
1
,
0
)
=
rdm
(
2
,
3
)
=
rdm
(
3
,
2
)
=
1
;
break
;
case
7
:
rdm
(
0
,
3
)
=
rdm
(
2
,
1
)
=
1
;
rdm
(
1
,
2
)
=
rdm
(
3
,
0
)
=
-
1
;
break
;
case
8
:
rdm
(
0
,
1
)
=
rdm
(
3
,
2
)
=
1
;
rdm
(
1
,
0
)
=
rdm
(
2
,
3
)
=
-
1
;
break
;
case
9
:
rdm
(
0
,
2
)
=
rdm
(
1
,
3
)
=
-
1
;
rdm
(
2
,
0
)
=
rdm
(
3
,
1
)
=
1
;
break
;
case
0
:
rdm
(
0
,
1
)
=
rdm
(
2
,
3
)
=
1
;
rdm
(
1
,
0
)
=
rdm
(
3
,
2
)
=
-
1
;
break
;
case
1
:
rdm
(
0
,
1
)
=
rdm
(
1
,
0
)
=
-
1
;
rdm
(
2
,
3
)
=
rdm
(
3
,
2
)
=
1
;
break
;
case
2
:
rdm
(
0
,
3
)
=
rdm
(
1
,
2
)
=
rdm
(
2
,
1
)
=
rdm
(
3
,
0
)
=
1
;
break
;
case
3
:
rdm
(
0
,
0
)
=
rdm
(
2
,
2
)
=
-
1
;
rdm
(
1
,
1
)
=
rdm
(
3
,
3
)
=
1
;
break
;
case
4
:
rdm
(
0
,
0
)
=
rdm
(
3
,
3
)
=
-
1
;
rdm
(
1
,
1
)
=
rdm
(
2
,
2
)
=
1
;
break
;
case
5
:
rdm
(
0
,
2
)
=
rdm
(
2
,
0
)
=
1
;
rdm
(
1
,
3
)
=
rdm
(
3
,
1
)
=
-
1
;
break
;
case
6
:
rdm
(
0
,
1
)
=
rdm
(
1
,
0
)
=
rdm
(
2
,
3
)
=
rdm
(
3
,
2
)
=
1
;
break
;
case
7
:
rdm
(
0
,
3
)
=
rdm
(
2
,
1
)
=
1
;
rdm
(
1
,
2
)
=
rdm
(
3
,
0
)
=
-
1
;
break
;
case
8
:
rdm
(
0
,
1
)
=
rdm
(
3
,
2
)
=
1
;
rdm
(
1
,
0
)
=
rdm
(
2
,
3
)
=
-
1
;
break
;
case
9
:
rdm
(
0
,
2
)
=
rdm
(
1
,
3
)
=
-
1
;
rdm
(
2
,
0
)
=
rdm
(
3
,
1
)
=
1
;
break
;
case
10
:
rdm
(
0
,
2
)
=
rdm
(
3
,
1
)
=
1
;
rdm
(
1
,
3
)
=
rdm
(
2
,
0
)
=
-
1
;
break
;
case
11
:
rdm
(
0
,
2
)
=
rdm
(
1
,
3
)
=
rdm
(
2
,
0
)
=
rdm
(
3
,
1
)
=
-
1
;
break
;
case
12
:
rdm
(
0
,
0
)
=
rdm
(
1
,
1
)
=
-
1
;
rdm
(
2
,
2
)
=
rdm
(
3
,
3
)
=
1
;
break
;
case
13
:
rdm
(
0
,
3
)
=
rdm
(
3
,
0
)
=
-
1
;
rdm
(
1
,
2
)
=
rdm
(
2
,
1
)
=
1
;
break
;
case
14
:
rdm
(
0
,
3
)
=
rdm
(
1
,
2
)
=
-
1
;
rdm
(
2
,
1
)
=
rdm
(
3
,
0
)
=
1
;
break
;
case
15
:
rdm
(
0
,
0
)
=
rdm
(
1
,
1
)
=
rdm
(
2
,
2
)
=
rdm
(
3
,
3
)
=
1
;
break
;
default:
throw
OpalException
(
"R
DM
::getRDM()"
,
default:
throw
OpalException
(
"R
ealDiracMatrix
::getRDM()"
,
"Index (i = "
+
std
::
to_string
(
i
)
+
" out of range: 0 <= i <= 15"
);
break
;
}
...
...
@@ -57,14 +58,14 @@ RDM::getRDM(short i) {
}
typename
R
DM
::
vector_t
R
DM
::
decompose
(
const
matrix_t
&
M
)
{
typename
R
ealDiracMatrix
::
vector_t
R
ealDiracMatrix
::
decompose
(
const
matrix_t
&
M
)
{
/*
* Formula (11) from paper:
* Geometrical method of decoupling
*/
if
(
M
.
size1
()
!=
4
&&
M
.
size1
()
!=
M
.
size2
())
{
throw
OpalException
(
"R
DM
::decompose()"
,
throw
OpalException
(
"R
ealDiracMatrix
::decompose()"
,
"Wrong matrix dimensions."
);
}
...
...
@@ -89,10 +90,10 @@ RDM::decompose(const matrix_t& M) {
}
typename
R
DM
::
matrix_t
R
DM
::
combine
(
const
vector_t
&
coeffs
)
{
typename
R
ealDiracMatrix
::
matrix_t
R
ealDiracMatrix
::
combine
(
const
vector_t
&
coeffs
)
{
if
(
coeffs
.
size
()
>
16
)
{
throw
OpalException
(
"R
DM
::combine()"
,
throw
OpalException
(
"R
ealDiracMatrix
::combine()"
,
"Wrong size of coefficient vector."
);
}
...
...
@@ -107,7 +108,9 @@ RDM::combine(const vector_t& coeffs) {
}
void
RDM
::
diagonalize
(
matrix_t
&
Ms
,
sparse_matrix_t
&
R
,
sparse_matrix_t
&
invR
)
{
#include <iomanip>
void
RealDiracMatrix
::
diagonalize
(
matrix_t
&
Ms
,
sparse_matrix_t
&
R
,
sparse_matrix_t
&
invR
)
{
// R and invR store the total transformation
R
=
boost
::
numeric
::
ublas
::
identity_matrix
<
double
>
(
4
);
...
...
@@ -127,36 +130,36 @@ void RDM::diagonalize(matrix_t& Ms, sparse_matrix_t& R, sparse_matrix_t& invR) {
};
// 1. Transformation with \gamma_{0}
P
(
0
)
=
mult
(
1
);
P
(
1
)
=
mult
(
2
);
P
(
2
)
=
mult
(
3
);
E
(
0
)
=
mult
(
4
);
E
(
1
)
=
mult
(
5
);
E
(
2
)
=
mult
(
6
);
P
(
0
)
=
mult
(
1
);
P
(
1
)
=
mult
(
2
);
P
(
2
)
=
mult
(
3
);
E
(
0
)
=
mult
(
4
);
E
(
1
)
=
mult
(
5
);
E
(
2
)
=
mult
(
6
);
B
(
0
)
=
-
mult
(
7
);
B
(
1
)
=
-
mult
(
8
);
B
(
2
)
=
-
mult
(
9
);
mr
=
boost
::
numeric
::
ublas
::
inner_prod
(
E
,
B
);
// formula (31), paper: Geometrical method of decoupling
mg
=
boost
::
numeric
::
ublas
::
inner_prod
(
B
,
P
);
// formula (31), paper: Geometrical method of decoupling
transform
(
Ms
,
short
(
0
),
0.5
*
std
::
atan2
(
mg
,
mr
),
R
,
invR
);
transform
(
Ms
,
0
,
0.5
*
std
::
atan2
(
mg
,
mr
),
R
,
invR
);
// 2. Transformation with \gamma_{7}
eps
=
-
mult
(
0
);
P
(
0
)
=
mult
(
1
);
P
(
1
)
=
mult
(
2
);
P
(
2
)
=
mult
(
3
);
E
(
0
)
=
mult
(
4
);
E
(
1
)
=
mult
(
5
);
E
(
2
)
=
mult
(
6
);
eps
=
-
mult
(
0
);
P
(
0
)
=
mult
(
1
);
P
(
1
)
=
mult
(
2
);
P
(
2
)
=
mult
(
3
);
E
(
0
)
=
mult
(
4
);
E
(
1
)
=
mult
(
5
);
E
(
2
)
=
mult
(
6
);
B
(
0
)
=
-
mult
(
7
);
B
(
1
)
=
-
mult
(
8
);
B
(
2
)
=
-
mult
(
9
);
b
=
eps
*
B
+
matt_boost
::
cross_prod
(
E
,
P
);
// formula (32), paper: Geometrical method of decoupling
transform
(
Ms
,
short
(
7
),
0.5
*
std
::
atan2
(
b
(
2
),
b
(
1
)),
R
,
invR
);
transform
(
Ms
,
7
,
0.5
*
std
::
atan2
(
b
(
2
),
b
(
1
)),
R
,
invR
);
// 3. Transformation with \gamma_{9}
eps
=
-
mult
(
0
);
P
(
0
)
=
mult
(
1
);
P
(
1
)
=
mult
(
2
);
P
(
2
)
=
mult
(
3
);
E
(
0
)
=
mult
(
4
);
E
(
1
)
=
mult
(
5
);
E
(
2
)
=
mult
(
6
);
eps
=
-
mult
(
0
);
P
(
0
)
=
mult
(
1
);
P
(
1
)
=
mult
(
2
);
P
(
2
)
=
mult
(
3
);
E
(
0
)
=
mult
(
4
);
E
(
1
)
=
mult
(
5
);
E
(
2
)
=
mult
(
6
);
B
(
0
)
=
-
mult
(
7
);
B
(
1
)
=
-
mult
(
8
);
B
(
2
)
=
-
mult
(
9
);
b
=
eps
*
B
+
matt_boost
::
cross_prod
(
E
,
P
);
transform
(
Ms
,
short
(
9
),
-
0.5
*
std
::
atan2
(
b
(
0
),
b
(
1
)),
R
,
invR
);
transform
(
Ms
,
9
,
-
0.5
*
std
::
atan2
(
b
(
0
),
b
(
1
)),
R
,
invR
);
// 4. Transformation with \gamma_{2}
eps
=
-
mult
(
0
);
P
(
0
)
=
mult
(
1
);
P
(
1
)
=
mult
(
2
);
P
(
2
)
=
mult
(
3
);
E
(
0
)
=
mult
(
4
);
E
(
1
)
=
mult
(
5
);
E
(
2
)
=
mult
(
6
);
eps
=
-
mult
(
0
);
P
(
0
)
=
mult
(
1
);
P
(
1
)
=
mult
(
2
);
P
(
2
)
=
mult
(
3
);
E
(
0
)
=
mult
(
4
);
E
(
1
)
=
mult
(
5
);
E
(
2
)
=
mult
(
6
);
B
(
0
)
=
-
mult
(
7
);
B
(
1
)
=
-
mult
(
8
);
B
(
2
)
=
-
mult
(
9
);
mr
=
boost
::
numeric
::
ublas
::
inner_prod
(
E
,
B
);
b
=
eps
*
B
+
matt_boost
::
cross_prod
(
E
,
P
);
...
...
@@ -165,14 +168,14 @@ void RDM::diagonalize(matrix_t& Ms, sparse_matrix_t& R, sparse_matrix_t& invR) {
// in rdm.c of Dr. Christian Baumgarten
if
(
std
::
fabs
(
mr
)
<
std
::
fabs
(
b
(
1
)))
{
transform
(
Ms
,
short
(
2
),
0.5
*
std
::
atanh
(
mr
/
b
(
1
)),
R
,
invR
);
transform
(
Ms
,
2
,
0.5
*
std
::
atanh
(
mr
/
b
(
1
)),
R
,
invR
);
}
else
{
transform
(
Ms
,
short
(
2
),
0.5
*
std
::
atanh
(
b
(
1
)
/
mr
),
R
,
invR
);
transform
(
Ms
,
2
,
0.5
*
std
::
atanh
(
b
(
1
)
/
mr
),
R
,
invR
);
}
eps
=
-
mult
(
0
);
P
(
0
)
=
mult
(
1
);
P
(
1
)
=
mult
(
2
);
P
(
2
)
=
mult
(
3
);
E
(
0
)
=
mult
(
4
);
E
(
1
)
=
mult
(
5
);
E
(
2
)
=
mult
(
6
);
eps
=
-
mult
(
0
);
P
(
0
)
=
mult
(
1
);
P
(
1
)
=
mult
(
2
);
P
(
2
)
=
mult
(
3
);
E
(
0
)
=
mult
(
4
);
E
(
1
)
=
mult
(
5
);
E
(
2
)
=
mult
(
6
);
B
(
0
)
=
-
mult
(
7
);
B
(
1
)
=
-
mult
(
8
);
B
(
2
)
=
-
mult
(
9
);
// formula (31), paper: Geometrical method of decoupling
...
...
@@ -184,17 +187,37 @@ void RDM::diagonalize(matrix_t& Ms, sparse_matrix_t& R, sparse_matrix_t& invR) {
double
E2
=
boost
::
numeric
::
ublas
::
inner_prod
(
E
,
E
);
// 5. Transform with \gamma_{0}
transform
(
Ms
,
short
(
0
),
0.25
*
std
::
atan2
(
mb
,
0.5
*
(
E2
-
P2
)),
R
,
invR
);
transform
(
Ms
,
0
,
0.25
*
std
::
atan2
(
mb
,
0.5
*
(
E2
-
P2
)),
R
,
invR
);
// 6. Transformation with \gamma_{8}
P
(
0
)
=
mult
(
1
);
P
(
2
)
=
mult
(
3
);
transform
(
Ms
,
short
(
8
),
-
0.5
*
std
::
atan2
(
P
(
2
),
P
(
0
)),
R
,
invR
);
transform
(
Ms
,
8
,
-
0.5
*
std
::
atan2
(
P
(
2
),
P
(
0
)),
R
,
invR
);
}
RealDiracMatrix
::
sparse_matrix_t
RealDiracMatrix
::
diagonalize
(
matrix_t
&