~maus-maintainers/maus/merge

1
2
3
4
5
6
7
8
9
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
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
/* This file is part of MAUS: http://   micewww.pp.rl.ac.uk:8080/projects/maus
 * 
 * MAUS 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.
 * 
 * MAUS is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with MAUS.  If not, see <http://   www.gnu.org/licenses/>.
 */

/* Author: Peter Lane
 */

#include "src/common_cpp/Optics/CovarianceMatrix.hh"

#include <cmath>
#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <sstream>

#include "CLHEP/Matrix/SymMatrix.h"
#include "CLHEP/Units/PhysicalConstants.h"
#include "TMatrixDSym.h"

#include "Utils/Exception.hh"
#include "Maths/Matrix.hh"
#include "Maths/SymmetricMatrix.hh"

namespace MAUS {

using MAUS::Matrix;
using ::CLHEP::c_light;

// ############################
//  Free Functions
// ############################

// ****************************
//  Conversion Functions
// ****************************

CovarianceMatrix rotate(const CovarianceMatrix& covariances,
                        const double angle) {
  double fcos = ::cos(angle);
  double fsin = ::sin(angle);
  double rotation_array[36] = {
    1.,    0.,    0.,    0.,    0.,    0.,
    0.,    1.,    0.,    0.,    0.,    0.,
    0.,    0.,    fcos,  0.,    fsin,  0.,
    0.,    0.,    0.,    fcos,  0.,    fsin,
    0.,    0.,    -fsin, 0.,    fcos,  0.,
    0.,    0.,    0.,    -fsin, 0.,    fcos
  };
  const Matrix<double> rotation(6, 6, rotation_array);
  const Matrix<double> rotation_transpose = transpose(rotation);

  // the orthogonal similarity transform of a symmetric matrix is also symmetric
  CovarianceMatrix rotated_covariances(
    rotation * covariances * rotation_transpose);

  return rotated_covariances;
}

// ############################
//  CovarianceMatrix
// ############################

// ****************************
//  Constructors
// ****************************

CovarianceMatrix::CovarianceMatrix() : SymmetricMatrix(6)
{ }

CovarianceMatrix::CovarianceMatrix(const CovarianceMatrix& original_instance)
    : SymmetricMatrix((SymmetricMatrix) original_instance)
{ }

CovarianceMatrix::CovarianceMatrix(const Matrix<double>& matrix)
    : SymmetricMatrix() {
  if (   (matrix.number_of_rows() < 6)
      || (matrix.number_of_columns() < 6)) {
    throw(Exception(Exception::recoverable,
                 "Attempted to construct with a Matrix<double> containing "
                 "fewer than six rows and/or columns",
                 "CovarianceMatrix::CovarianceMatrix(Matrix<double>)"));
  }
  build_matrix(6);
  for (size_t row = 1; row <= 6; ++row) {
    for (size_t column = 1; column <= row; ++column) {
      set(row, column, matrix(row, column));
    }
  }
}

CovarianceMatrix::CovarianceMatrix(const SymmetricMatrix& symmetric_matrix)
    : SymmetricMatrix() {
  if (   (symmetric_matrix.number_of_rows() < 6)
      || (symmetric_matrix.number_of_columns() < 6)) {
    throw(Exception(Exception::recoverable,
          "Attempted to construct with a SymmetricMatrix containing "
          "fewer than six rows/columns",
          "CovarianceMatrix::CovarianceMatrix(SymmetricMatrix<double>)"));
  }
  build_matrix(6);
  for (size_t row = 1; row <= 6; ++row) {
    for (size_t column = 1; column <= row; ++column) {
      set(row, column, symmetric_matrix(row, column));
    }
  }
}

CovarianceMatrix::CovarianceMatrix(const ::CLHEP::HepSymMatrix& hep_sym_matrix)
    : SymmetricMatrix() {
  if (   (hep_sym_matrix.num_row() < 6)
      || (hep_sym_matrix.num_col() < 6)) {
    throw(Exception(Exception::recoverable,
          "Attempted to construct with a HepSymMatrix containing fewer "
          "than six rows/columns",
          "CovarianceMatrix::CovarianceMatrix(CLHEP::HepSymMatrix<double>)"));
  }
  build_matrix(6);
  double element;
  for (size_t row = 1; row <= 6; ++row) {
    for (size_t column = 1; column <= row; ++column) {
      element = hep_sym_matrix(row, column);
      Matrix<double>::operator()(row, column) = element;
      Matrix<double>::operator()(column, row) = element;
    }
  }
}

CovarianceMatrix::CovarianceMatrix(const TMatrixDSym& root_sym_matrix)
    : SymmetricMatrix(root_sym_matrix) {
  if (   (root_sym_matrix.GetNrows() < 6)
      || (root_sym_matrix.GetNcols() < 6)) {
    throw(Exception(Exception::recoverable,
          "Attempted to construct with a TMatrixDSym containing fewer "
          "than six rows/columns",
          "CovarianceMatrix::CovarianceMatrix(TMatrixDSym<double>)"));
  }
  const double * data = root_sym_matrix.GetMatrixArray();
  build_matrix(6, data);
}

CovarianceMatrix::CovarianceMatrix(double const * const array_matrix)
    : SymmetricMatrix(6, array_matrix)
{ }

/* xboa uses an equation for calculating canonical angular momentum that is
 * different from the standard. This in turn affects Ltwiddle_t. Furthermore,
 * xboa uses a different definition of kappa and sigma_y_Px than G. Penn's
 * paper (this might be an error). It also appears that xboa uses total
 * momentum instead of longitudinal momentum in calculations of the covariances.
 */
const CovarianceMatrix CovarianceMatrix::CreateFromPennParameters(
    double mass,
    double momentum,
    double charge,
    double Bz,
    double Ltwiddle_t,
    double emittance_t,
    double beta_t,
    double alpha_t,
    double emittance_l,
    double beta_l,
    double alpha_l,
    double dispersion_x,
    double dispersion_prime_x,
    double dispersion_y,
    double dispersion_prime_y) {

  // *** calculate some intermediate values ***
  double energy     = ::sqrt(momentum * momentum + mass * mass);
  // this differs from G. Penn's paper
  double kappa      = c_light * Bz / (2. * momentum);
  double gamma_t    = (1 + alpha_t*alpha_t
                       + (beta_t*kappa - Ltwiddle_t)*(beta_t*kappa-Ltwiddle_t))
                    / beta_t;
  double gamma_l    = (1+alpha_l*alpha_l+beta_l*beta_l*kappa*kappa)/beta_l;

  // *** calculate the lower triangle covariances <A B> = sigma_A_B ***
  double sigma_t_t  =  emittance_l * mass * beta_l / momentum;

  double sigma_E_t  = -emittance_l * mass * alpha_l;
  double sigma_E_E  =  emittance_l * mass * gamma_l * momentum;

  double sigma_x_t  =  0.;
  double sigma_x_E  =  dispersion_x * sigma_E_E / energy;
  double sigma_x_x  =  emittance_t * beta_t * mass / momentum;

  double sigma_Px_t =  0.;
  double sigma_Px_E =  dispersion_prime_x * sigma_E_E / energy;
  double sigma_Px_x = -mass * emittance_t * alpha_t;
  double sigma_Px_Px=  mass * momentum * emittance_t * gamma_t;

  double sigma_y_t  =  0.;
  double sigma_y_E  = dispersion_y * sigma_E_E / energy;
  double sigma_y_x  =   0.;
  // this differes from G. Penn's paper
  double sigma_y_Px =  mass * emittance_t
                    * (beta_t*kappa - Ltwiddle_t) * charge;
  double sigma_y_y  =  emittance_t * beta_t * mass / momentum;

  double sigma_Py_t =  0.;
  double sigma_Py_E =  dispersion_prime_y * sigma_E_E / energy;
  double sigma_Py_x = -sigma_y_Px;
  double sigma_Py_Px=  0.;
  double sigma_Py_y =  sigma_Px_x;
  double sigma_Py_Py=  sigma_Px_Px;

  double covariances[36] = {
    sigma_t_t,  0.,         0.,         0.,           0.,         0.,
    sigma_E_t,  sigma_E_E,  0.,         0.,           0.,         0.,
    sigma_x_t,  sigma_x_E,  sigma_x_x,  0.,           0.,         0.,
    sigma_Px_t, sigma_Px_E, sigma_Px_x, sigma_Px_Px,  0.,         0.,
    sigma_y_t,  sigma_y_E,  sigma_y_x,  sigma_y_Px,   sigma_y_y,  0.,
    sigma_Py_t, sigma_Py_E, sigma_Py_x, sigma_Py_Px,  sigma_Py_y, sigma_Py_Py
  };

  return CovarianceMatrix(covariances);
}

const CovarianceMatrix CovarianceMatrix::CreateFromTwissParameters(
      double mass,
      double momentum,
      double emittance_x,
      double beta_x,
      double alpha_x,
      double emittance_y,
      double beta_y,
      double alpha_y,
      double emittance_l,
      double beta_l,
      double alpha_l,
      double dispersion_x,
      double dispersion_prime_x,
      double dispersion_y,
      double dispersion_prime_y) {
  // *** calculate some intermediate values ***
  double energy  = ::sqrt(momentum * momentum + mass * mass);
  double gamma_x = (1+alpha_x*alpha_x)/beta_x;
  double gamma_y = (1+alpha_y*alpha_y)/beta_y;
  double gamma_l = (1+alpha_l*alpha_l)/beta_l;

  // *** calculate the lower triangle covariances <A B> = sigma_A_B ***
  double sigma_t_t  = emittance_l * mass * beta_l / momentum;

  double sigma_E_t  = -emittance_l * mass * alpha_l;
  double sigma_E_E  = emittance_l * mass * gamma_l * momentum;

  double sigma_x_t  = 0.;
  double sigma_x_E  = dispersion_x * sigma_E_E / energy;
  // FIXME(plane1@hawk.iit.edu) Shouldn't this be just beta_x emittance_x?
  double sigma_x_x  = emittance_x * mass * beta_x / momentum;

  double sigma_Px_t = 0.;
  double sigma_Px_E = dispersion_prime_x * sigma_E_E / energy;
  // FIXME(plane1@hawk.iit.edu) Shouldn't this be just -alpha_x emittance_x ?
  double sigma_Px_x = -emittance_x * mass * alpha_x;
  // FIXME(plane1@hawk.iit.edu) Shouldn't this be just gamma_x emittance_x ?
  double sigma_Px_Px= emittance_x * mass * momentum * gamma_x;

  double sigma_y_t  = 0.;
  double sigma_y_E  = dispersion_y * sigma_E_E / energy;
  double sigma_y_x  = 0.;
  double sigma_y_Px = 0.;
  // FIXME(plane1@hawk.iit.edu) Shouldn't this be just beta_y emittance_y?
  double sigma_y_y  = emittance_y * mass * beta_y / momentum;

  double sigma_Py_t = 0.;
  double sigma_Py_E = dispersion_prime_y * sigma_E_E / energy;
  double sigma_Py_x = 0.;
  double sigma_Py_Px= 0.;
  // FIXME(plane1@hawk.iit.edu) Shouldn't this be just -alpha_y emittance_y ?
  double sigma_Py_y = -emittance_y * mass * alpha_y;
  // FIXME(plane1@hawk.iit.edu) Shouldn't this be just gamma_y emittance_y ?
  double sigma_Py_Py= emittance_y * mass * momentum * gamma_y;

  double covariances[36] = {
    sigma_t_t,  0.,         0.,         0.,           0.,         0.,
    sigma_E_t,  sigma_E_E,  0.,         0.,           0.,         0.,
    sigma_x_t,  sigma_x_E,  sigma_x_x,  0.,           0.,         0.,
    sigma_Px_t, sigma_Px_E, sigma_Px_x, sigma_Px_Px,  0.,         0.,
    sigma_y_t,  sigma_y_E,  sigma_y_x,  sigma_y_Px,   sigma_y_y,  0.,
    sigma_Py_t, sigma_Py_E, sigma_Py_x, sigma_Py_Px,  sigma_Py_y, sigma_Py_Py
  };

  return CovarianceMatrix(covariances);
}

bool CovarianceMatrix::IsPositiveDefinite() const {
  size_t min_row = 1;
  size_t min_column = 1;
  size_t rows = size();
  for (size_t length = 1; length <= rows; length++) {
    // Sylvester Criterion - all the leading principle minors must be positive
    if (determinant(submatrix(min_row, length, min_column, length)) <= 0) {
      return false;
    }
  }

  return true;
}
}  // namespace MAUS