Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Classes | Public Types | Public Member Functions | Static Public Member Functions | Protected Member Functions | Protected Attributes | Static Private Member Functions | Private Attributes | List of all members
JFIT::JMatrixNZ Class Reference

Determination of the co-variance matrix of hits for a track along z-axis (JFIT::JLine1Z). More...

#include <JMatrixNZ.hh>

Inheritance diagram for JFIT::JMatrixNZ:
JMATH::JMatrixNS JMATH::JMatrixND JMATH::JMatrixND_t JMATH::JMath< JMatrixND > JLANG::JEquals< JFirst_t, JSecond_t > JLANG::JSingleton< T >

Classes

struct  variance
 Auxiliary data structure for co-variance calculation. More...
 

Public Types

typedef std::vector< std::pair
< size_t, size_t > > 
JPermutationMatrix
 Type definition of permutation matrix. More...
 
typedef T data_type
 

Public Member Functions

 JMatrixNZ ()
 Default contructor. More...
 
template<class T >
 JMatrixNZ (const JVector3D &pos, T __begin, T __end, const double alpha, const double sigma)
 Constructor. More...
 
template<class T >
void set (const JVector3D &pos, T __begin, T __end, const double alpha, const double sigma)
 Set co-variance matrix. More...
 
void invert ()
 Invert matrix according LDU decomposition. More...
 
void update (const size_t k, const double value)
 Update inverted matrix at given diagonal element. More...
 
void resize (const size_t size)
 Resize matrix. More...
 
JMatrixNDreset ()
 Set matrix to the null matrix. More...
 
JMatrixNDsetIdentity ()
 Set to identity matrix. More...
 
JMatrixNDnegate ()
 Negate matrix. More...
 
JMatrixNDadd (const JMatrixND &A)
 Matrix addition. More...
 
JMatrixNDsub (const JMatrixND &A)
 Matrix subtraction. More...
 
JMatrixNDmul (const double factor)
 Scale matrix. More...
 
JMatrixNDmul (const JMatrixND &A, const JMatrixND &B)
 Matrix multiplication. More...
 
JMatrixNDmul (const JNullType &object)
 Multiply with object. More...
 
JMatrixNDdiv (const double factor)
 Scale matrix. More...
 
bool equals (const JMatrixND &A, const double eps=std::numeric_limits< double >::min()) const
 Equality. More...
 
bool isIdentity (const double eps=std::numeric_limits< double >::min()) const
 Test identity. More...
 
bool isSymmetric (const double eps=std::numeric_limits< double >::min()) const
 Test symmetry. More...
 
double getDot (const JVectorND &v) const
 Get dot product. More...
 
void clear ()
 Clear memory. More...
 
void set (const JMatrixND_t &A)
 Set matrix. More...
 
void swap (JMatrixND_t &A)
 Swap matrices. More...
 
JMatrixND_ttranspose ()
 Transpose. More...
 
size_t size () const
 Get dimension of matrix. More...
 
size_t capacity () const
 Get capacity of dimension. More...
 
bool empty () const
 Check emptyness. More...
 
const double * data () const
 Get pointer to data. More...
 
double * data ()
 Get pointer to data. More...
 
const double * operator[] (size_t row) const
 Get row data. More...
 
double * operator[] (size_t row)
 Get row data. More...
 
double operator() (const size_t row, const size_t col) const
 Get matrix element. More...
 
double & operator() (const size_t row, const size_t col)
 Get matrix element. More...
 
double at (size_t row, size_t col) const
 Get matrix element. More...
 
double & at (size_t row, size_t col)
 Get matrix element. More...
 

Static Public Member Functions

static data_typegetInstance ()
 Get unique instance of template class. More...
 

Protected Member Functions

void rswap (size_t r1, size_t r2)
 
void cswap (size_t c1, size_t c2)
 Swap columns. More...
 

Protected Attributes

double * __p
 pointer to data More...
 
size_t __n
 dimension of matrix More...
 
size_t __m
 capacity of matrix More...
 

Static Private Member Functions

static double getDot (const variance &first, const variance &second)
 Get dot product. More...
 

Private Attributes

std::vector< variancebuffer
 

Detailed Description

Determination of the co-variance matrix of hits for a track along z-axis (JFIT::JLine1Z).

In this, the given angular and time resolution are taken into account.

Definition at line 28 of file JMatrixNZ.hh.

Member Typedef Documentation

typedef std::vector< std::pair<size_t,size_t> > JMATH::JMatrixNS::JPermutationMatrix
inherited

Type definition of permutation matrix.

Definition at line 33 of file JMatrixNS.hh.

template<class T>
typedef T JLANG::JSingleton< T >::data_type
inherited

Definition at line 20 of file JSingleton.hh.

Constructor & Destructor Documentation

JFIT::JMatrixNZ::JMatrixNZ ( )
inline

Default contructor.

Definition at line 39 of file JMatrixNZ.hh.

39  :
40  JMatrixNS()
41  {}
JMatrixNS()
Default constructor.
Definition: JMatrixNS.hh:39
template<class T >
JFIT::JMatrixNZ::JMatrixNZ ( const JVector3D pos,
T  __begin,
T  __end,
const double  alpha,
const double  sigma 
)
inline

Constructor.

The template argument T refers to an iterator of a data structure which should have the following member methods:

  • double getX(); // [m]
  • double getY(); // [m]
  • double getZ(); // [m]
Parameters
posreference position [m]
__beginbegin of data
__endend of data
alphaangular resolution [deg]
sigmatime resolution [ns]

Definition at line 59 of file JMatrixNZ.hh.

63  :
64  JMatrixNS()
65  {
66  set(pos, __begin, __end, alpha, sigma);
67  }
void set(const JVector3D &pos, T __begin, T __end, const double alpha, const double sigma)
Set co-variance matrix.
Definition: JMatrixNZ.hh:85
JMatrixNS()
Default constructor.
Definition: JMatrixNS.hh:39

Member Function Documentation

template<class T >
void JFIT::JMatrixNZ::set ( const JVector3D pos,
T  __begin,
T  __end,
const double  alpha,
const double  sigma 
)
inline

Set co-variance matrix.

The template argument T refers to an iterator of a data structure which should have the following member methods:

  • double getX(); // [m]
  • double getY(); // [m]
  • double getZ(); // [m]
Parameters
posreference position [m]
__beginbegin of data
__endend of data
alphaangular resolution [deg]
sigmatime resolution [ns]

Definition at line 85 of file JMatrixNZ.hh.

90  {
91  using namespace std;
92  using namespace JTOOLS;
93 
94 
95  const int N = distance(__begin, __end);
96 
97  this ->resize(N);
98  buffer.resize(N);
99 
100  const double ta = alpha * PI / 180.0;
101  const double ct = cos(ta);
102  const double st = sin(ta);
103 
104 
105  // angular resolution
106 
107  for (T i = __begin; i != __end; ++i) {
108 
109  const double dx = i->getX() - pos.getX();
110  const double dy = i->getY() - pos.getY();
111  const double dz = i->getZ() - pos.getZ();
112 
113  const double R = sqrt(dx*dx + dy*dy);
114 
115  double x = ta * getKappaC() * getInverseSpeedOfLight();
116  double y = ta * getKappaC() * getInverseSpeedOfLight();
117  double v = ta * getInverseSpeedOfLight();
118  double w = ta * getInverseSpeedOfLight();
119 
120  if (R != 0.0) {
121  x *= dx / R;
122  y *= dy / R;
123  }
124 
125  x *= (dz*ct - dx*st);
126  y *= (dz*ct - dy*st);
127  v *= -(dx*ct + dz*st);
128  w *= -(dy*ct + dz*st);
129 
130  buffer[distance(__begin,i)] = variance(x, y, v, w);
131  }
132 
133  for (int i = 0; i != N; ++i) {
134 
135  for (int j = 0; j != i; ++j) {
136  (*this)(i, j) = getDot(buffer[i], buffer[j]);
137  (*this)(j, i) = (*this)(i, j);
138  }
139 
140  (*this)(i, i) = getDot(buffer[i], buffer[i]) + sigma * sigma;
141  }
142  }
data_type w[N+1][M+1]
Definition: JPolint.hh:708
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
static double getDot(const variance &first, const variance &second)
Get dot product.
Definition: JMatrixNZ.hh:196
static const double PI
Constants.
Definition: JConstants.hh:20
const double getInverseSpeedOfLight()
Get inverse speed of light.
Definition: JConstants.hh:100
do set_variable OUTPUT_DIRECTORY $WORKDIR T
std::vector< variance > buffer
Definition: JMatrixNZ.hh:186
then usage $script[distance] fi case set_variable R
Definition: JDrawLED.sh:40
int j
Definition: JPolint.hh:634
data_type v[N+1][M+1]
Definition: JPolint.hh:707
double getKappaC()
Get average kappa of Cherenkov light in water.
Definition: JConstants.hh:166
then usage $script[input file[working directory[option]]] nWhere option can be N
Definition: JMuonPostfit.sh:37
static double JFIT::JMatrixNZ::getDot ( const variance first,
const variance second 
)
inlinestaticprivate

Get dot product.

Parameters
firstfirst variance
secondsecond variance
Returns
dot product

Definition at line 196 of file JMatrixNZ.hh.

197  {
198  return (first.x * second.x +
199  first.y * second.y +
200  first.v * second.v +
201  first.w * second.w);
202 
203  }
then echo The file $DIR KM3NeT_00000001_00000000 root already please rename or remove it first
void JMATH::JMatrixNS::invert ( )
inlineinherited

Invert matrix according LDU decomposition.

Definition at line 80 of file JMatrixNS.hh.

81  {
82  using namespace std;
83 
84  size_t i, row, col, root;
85 
86  double* p0;
87  double* p1;
88  double* p2;
89  double* p3;
90 
91  const double* c0;
92  const double* c1;
93  const double* c2;
94  const double* c3;
95 
96  const double* a0;
97  const double* a1;
98  const double* a2;
99  const double* a3;
100 
101  double v0;
102  double v1;
103  double v2;
104  double v3;
105 
106  double w0;
107  double w1;
108  double w2;
109  double w3;
110 
112 
113 
114  // LDU decomposition
115 
116  for (root = 0; root != this->size(); ++root) {
117 
118  double value = (*this)(root, root);
119  size_t pivot = root;
120 
121  for (row = root; ++row != this->size(); ) {
122  if (fabs((*this)(row, root)) > fabs(value)) {
123  value = (*this)(row, root);
124  pivot = row;
125  }
126  }
127 
128  if (value == 0.0) {
129  THROW(JDivisionByZero, "LDU decomposition zero pivot at " << root);
130  }
131 
132  if (pivot != root) {
133 
134  this->rswap(root, pivot);
135 
136  P.push_back(make_pair(root, pivot));
137  }
138 
139  for (row = root + 1; row + 4 <= this->size(); row += 4) { // process rows by four
140 
141  p0 = (*this)[row + 0] + root;
142  p1 = (*this)[row + 1] + root;
143  p2 = (*this)[row + 2] + root;
144  p3 = (*this)[row + 3] + root;
145 
146  v0 = (*p0++ /= value);
147  v1 = (*p1++ /= value);
148  v2 = (*p2++ /= value);
149  v3 = (*p3++ /= value);
150 
151  c0 = (*this)[root] + root + 1;
152 
153  for (i = root + 1; i + 4 <= this->size(); i += 4, p0 += 4, p1 += 4, p2 += 4, p3 += 4, c0 += 4) {
154  p0[0] -= v0 * c0[0]; p0[1] -= v0 * c0[1]; p0[2] -= v0 * c0[2]; p0[3] -= v0 * c0[3];
155  p1[0] -= v1 * c0[0]; p1[1] -= v1 * c0[1]; p1[2] -= v1 * c0[2]; p1[3] -= v1 * c0[3];
156  p2[0] -= v2 * c0[0]; p2[1] -= v2 * c0[1]; p2[2] -= v2 * c0[2]; p2[3] -= v2 * c0[3];
157  p3[0] -= v3 * c0[0]; p3[1] -= v3 * c0[1]; p3[2] -= v3 * c0[2]; p3[3] -= v3 * c0[3];
158  }
159 
160  for ( ; i != this->size(); ++i, ++p0, ++p1, ++p2, ++p3, ++c0) {
161  *p0 -= v0 * (*c0);
162  *p1 -= v1 * (*c0);
163  *p2 -= v2 * (*c0);
164  *p3 -= v3 * (*c0);
165  }
166  }
167 
168  for ( ; row != this->size(); ++row) { // remaining rows
169 
170  p0 = (*this)[row + 0] + root;
171 
172  v0 = (*p0++ /= value);
173 
174  c0 = (*this)[root] + root + 1;
175 
176  for (i = root + 1; i + 4 <= this->size(); i += 4, p0 += 4, c0 += 4) {
177  *p0 -= v0 * c0[0]; p0[1] -= v0 * c0[1]; p0[2] -= v0 * c0[2]; p0[3] -= v0 * c0[3];
178  }
179 
180  for ( ; i != this->size(); ++i, ++p0, ++c0) {
181  *p0 -= v0 * (*c0);
182  }
183  }
184  }
185 
186 
187  // invert D
188 
189  for (size_t i = 0; i != this->size(); ++i) {
190 
191  double& pivot = (*this)(i,i);
192 
193  if (pivot == 0.0) {
194  THROW(JDivisionByZero, "Invert D zero pivot at " << i);
195  }
196 
197  pivot = 1.0 / pivot;
198  }
199 
200 
201  // invert U
202 
204 
205  A.set(*this);
206  A.transpose();
207 
208  for (row = 0; row + 4 <= this->size(); row += 4) { // process rows by four
209 
210  p0 = (*this)[row + 0] + row + 0;
211  p1 = (*this)[row + 1] + row + 1;
212  p2 = (*this)[row + 2] + row + 2;
213  p3 = (*this)[row + 3] + row + 3;
214 
215  v0 = *(p0++);
216  v1 = *(p1++);
217  v2 = *(p2++);
218  v3 = *(p3++);
219 
220  for (col = row + 1; col + 4 <= this->size(); ++col, ++p0, ++p1, ++p2, ++p3) {
221 
222  w0 = (*p0) * -v0;
223  w1 = (*p1) * -v1;
224  w2 = (*p2) * -v2;
225  w3 = (*p3) * -v3;
226 
227  c0 = (*this)[row + 0] + row + 0 + 1;
228  c1 = (*this)[row + 1] + row + 1 + 1;
229  c2 = (*this)[row + 2] + row + 2 + 1;
230  c3 = (*this)[row + 3] + row + 3 + 1;
231 
232  a0 = A[col + 0] + row + 0 + 1;
233  a1 = A[col + 1] + row + 1 + 1;
234  a2 = A[col + 2] + row + 2 + 1;
235  a3 = A[col + 3] + row + 3 + 1;
236 
237  for ( ; c0 + 4 < p0; c0 += 4, c1 += 4, c2 += 4, c3 += 4, a0 += 4, a1 += 4, a2 += 4, a3 += 4) {
238  w0 -= c0[0]*a0[0] + c0[1]*a0[1] + c0[2]*a0[2] + c0[3]*a0[3];
239  w1 -= c1[0]*a1[0] + c1[1]*a1[1] + c1[2]*a1[2] + c1[3]*a1[3];
240  w2 -= c2[0]*a2[0] + c2[1]*a2[1] + c2[2]*a2[2] + c2[3]*a2[3];
241  w3 -= c3[0]*a3[0] + c3[1]*a3[1] + c3[2]*a3[2] + c3[3]*a3[3];
242  }
243 
244  for ( ; c0 != p0; ++c0, ++c1, ++c2, ++c3, ++a0, ++a1, ++a2, ++a3) {
245  w0 -= (*c0) * (*a0);
246  w1 -= (*c1) * (*a1);
247  w2 -= (*c2) * (*a2);
248  w3 -= (*c3) * (*a3);
249  }
250 
251  *p0 = w0 * (*this)(col + 0, col + 0);
252  *p1 = w1 * (*this)(col + 1, col + 1);
253  *p2 = w2 * (*this)(col + 2, col + 2);
254  *p3 = w3 * (*this)(col + 3, col + 3);
255  }
256 
257  // last 3 columns
258 
259  w0 = (*p0) * -v0;
260  w1 = (*p1) * -v1;
261  w2 = (*p2) * -v2;
262 
263  c0 = (*this)[row + 0] + row + 0 + 1;
264  c1 = (*this)[row + 1] + row + 1 + 1;
265  c2 = (*this)[row + 2] + row + 2 + 1;
266 
267  a0 = A[col + 0] + row + 0 + 1;
268  a1 = A[col + 1] + row + 1 + 1;
269  a2 = A[col + 2] + row + 2 + 1;
270 
271  for ( ; c0 + 4 < p0; c0 += 4, c1 += 4, c2 += 4, c3 += 4, a0 += 4, a1 += 4, a2 += 4, a3 += 4) {
272  w0 -= c0[0]*a0[0] + c0[1]*a0[1] + c0[2]*a0[2] + c0[3]*a0[3];
273  w1 -= c1[0]*a1[0] + c1[1]*a1[1] + c1[2]*a1[2] + c1[3]*a1[3];
274  w2 -= c2[0]*a2[0] + c2[1]*a2[1] + c2[2]*a2[2] + c2[3]*a2[3];
275  }
276 
277  for ( ; c0 != p0; ++c0, ++c1, ++c2, ++a0, ++a1, ++a2) {
278  w0 -= (*c0) * (*a0);
279  w1 -= (*c1) * (*a1);
280  w2 -= (*c2) * (*a2);
281  }
282 
283  *p0 = w0 * (*this)(col + 0, col + 0);
284  *p1 = w1 * (*this)(col + 1, col + 1);
285  *p2 = w2 * (*this)(col + 2, col + 2);
286 
287  ++col; ++p0; ++p1;
288 
289  // last 2 columns
290 
291  w0 = (*p0) * -v0;
292  w1 = (*p1) * -v1;
293 
294  c0 = (*this)[row + 0] + row + 0 + 1;
295  c1 = (*this)[row + 1] + row + 1 + 1;
296 
297  a0 = A[col + 0] + row + 0 + 1;
298  a1 = A[col + 1] + row + 1 + 1;
299 
300  for ( ; c0 + 4 < p0; c0 += 4, c1 += 4, c2 += 4, c3 += 4, a0 += 4, a1 += 4, a2 += 4, a3 += 4) {
301  w0 -= c0[0]*a0[0] + c0[1]*a0[1] + c0[2]*a0[2] + c0[3]*a0[3];
302  w1 -= c1[0]*a1[0] + c1[1]*a1[1] + c1[2]*a1[2] + c1[3]*a1[3];
303  }
304 
305  for ( ; c0 != p0; ++c0, ++c1, ++a0, ++a1) {
306  w0 -= (*c0) * (*a0);
307  w1 -= (*c1) * (*a1);
308  }
309 
310  *p0 = w0 * (*this)(col + 0, col + 0);
311  *p1 = w1 * (*this)(col + 1, col + 1);
312 
313  ++col; ++p0;
314 
315  // last column
316 
317  w0 = (*p0) * -v0;
318 
319  c0 = (*this)[row + 0] + row + 0 + 1;
320 
321  a0 = A[col + 0] + row + 0 + 1;
322 
323  for ( ; c0 + 4 < p0; c0 += 4, c1 += 4, c2 += 4, c3 += 4, a0 += 4, a1 += 4, a2 += 4, a3 += 4) {
324  w0 -= c0[0]*a0[0] + c0[1]*a0[1] + c0[2]*a0[2] + c0[3]*a0[3];
325  }
326 
327  for ( ; c0 != p0; ++c0, ++a0) {
328  w0 -= (*c0) * (*a0);
329  }
330 
331  *p0 = w0 * (*this)(col + 0, col + 0);
332  }
333 
334  for ( ; row != this->size(); ++row) { // remaining rows
335 
336  p0 = (*this)[row + 0] + row + 0;
337 
338  v0 = *(p0++);
339 
340  for (col = row + 1; col != this->size(); ++col, ++p0) {
341 
342  w0 = *p0;
343 
344  w0 *= -v0;
345 
346  c0 = (*this)[row + 0] + row + 1;
347 
348  a0 = A[col + 0] + row + 0 + 1;
349 
350  for ( ; c0 != p0; ++c0, ++a0) {
351  w0 -= (*c0) * (*a0);
352  }
353 
354  *p0 = w0 * (*this)(col + 0, col + 0);
355  }
356  }
357 
358 
359  // solve A^-1 x L = U^-1
360 
361  double* buffer = A.data();
362 
363  for (col = this->size() - 1; col != 0; --col) {
364 
365  for (i = col; i < this->size(); ++i) {
366 
367  buffer[i] = (*this)(i, col - 1);
368 
369  (*this)(i, col - 1) = 0.0;
370  }
371 
372  for (row = 0; row + 4 < this->size(); row += 4) {
373 
374  p0 = (*this)[row + 0] + col - 1;
375  p1 = (*this)[row + 1] + col - 1;
376  p2 = (*this)[row + 2] + col - 1;
377  p3 = (*this)[row + 3] + col - 1;
378 
379  c0 = p0 + 1;
380  c1 = p1 + 1;
381  c2 = p2 + 1;
382  c3 = p3 + 1;
383 
384  a0 = buffer + col;
385 
386  w0 = 0.0;
387  w1 = 0.0;
388  w2 = 0.0;
389  w3 = 0.0;
390 
391  for (i = col; i + 4 <= this->size(); i += 4, c0 += 4, c1 += 4, c2 += 4, c3 += 4, a0 += 4) {
392  w0 += c0[0]*a0[0] + c0[1]*a0[1] + c0[2]*a0[2] + c0[3]*a0[3];
393  w1 += c1[0]*a0[0] + c1[1]*a0[1] + c1[2]*a0[2] + c1[3]*a0[3];
394  w2 += c2[0]*a0[0] + c2[1]*a0[1] + c2[2]*a0[2] + c2[3]*a0[3];
395  w3 += c3[0]*a0[0] + c3[1]*a0[1] + c3[2]*a0[2] + c3[3]*a0[3];
396  }
397 
398  for ( ; i != this->size(); ++i, ++c0, ++c1, ++c2, ++c3, ++a0) {
399  w0 += (*c0) * (*a0);
400  w1 += (*c1) * (*a0);
401  w2 += (*c2) * (*a0);
402  w3 += (*c3) * (*a0);
403  }
404 
405  *p0 -= w0;
406  *p1 -= w1;
407  *p2 -= w2;
408  *p3 -= w3;
409  }
410 
411  for ( ; row < this->size(); ++row) {
412 
413  p0 = (*this)[row] + col - 1;
414 
415  c0 = p0 + 1;
416 
417  a0 = buffer + col;
418 
419  w0 = 0.0;
420 
421  for (i = col; i + 4 <= this->size(); i += 4, c0 += 4, c1 += 4, c2 += 4, c3 += 4, a0 += 4) {
422  w0 += c0[0]*a0[0] + c0[1]*a0[1] + c0[2]*a0[2] + c0[3]*a0[3];
423  }
424 
425  for ( ; i != this->size(); ++i, ++c0, ++a0) {
426  w0 += (*c0) * (*a0);
427  }
428 
429  *p0 -= w0;
430  }
431  }
432 
433 
434  for (JPermutationMatrix::reverse_iterator p = P.rbegin(); p != P.rend(); ++p) {
435  cswap(p->first, p->second);
436  }
437  }
TPaveText * p1
then JPlot1D f $WORKDIR postfit[prefit\] root
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:670
void rswap(size_t r1, size_t r2)
Definition: JMatrixND.hh:792
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:154
TCanvas * c1
Global variables to handle mouse events.
std::vector< std::pair< size_t, size_t > > JPermutationMatrix
Type definition of permutation matrix.
Definition: JMatrixNS.hh:33
static data_type & getInstance()
Get unique instance of template class.
Definition: JSingleton.hh:27
void cswap(size_t c1, size_t c2)
Swap columns.
Definition: JMatrixND.hh:808
Basic NxN matrix.
Definition: JMatrixND.hh:35
source $JPP_DIR setenv csh $JPP_DIR eval JShellParser o a A
then $DIR JPlotNPE PDG P
Definition: JPlotNPE-PDG.sh:60
void JMATH::JMatrixNS::update ( const size_t  k,
const double  value 
)
inlineinherited

Update inverted matrix at given diagonal element.

If A is the original matrix and A' is such that for each (i,j) != (k,k):

      A'[i][j] = A[i][j]
      A'[k][k] = A[k][k] + value

then JMatrixNS::update(k, value) will modify the already inverted matrix of A such that it will be the inverse of A'.

See also arXiv.

Parameters
kindex of diagonal element
valuevalue to add

Definition at line 456 of file JMatrixNS.hh.

457  {
458  if (k >= this->size()) {
459  THROW(JIndexOutOfRange, "JMatrixNS::update(): index out of range " << k);
460  }
461 
462  size_t i, row;
463  double* col;
464 
465  const double del = (*this)(k,k);
466  const double din = 1.0 / (1.0 + value * del);
467 
468  double* root = (*this)[k];
469 
470  for (row = 0; row != this->size(); ++row) {
471 
472  if (row != k) {
473 
474  const double val = (*this)(row, k);
475 
476  for (i = 0, col = (*this)[row]; i != this->size(); ++i, ++col) {
477  if (i != k) {
478  (*col) -= val * root[i] * value * din;
479  }
480  }
481 
482  (*this)(row, k) *= din;
483  }
484  }
485 
486  for (i = 0, col = root; i != this->size(); ++i, ++col) {
487  if (i != k) {
488  (*col) *= din;
489  }
490  }
491 
492  root[k] = del * din;
493  }
*fatal Wrong number of arguments esac JCookie sh JRuns D $DETECTOR d sort n k
Definition: JRunrange.sh:16
then JPlot1D f $WORKDIR postfit[prefit\] root
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:670
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:154
void JMATH::JMatrixND::resize ( const size_t  size)
inlineinherited

Resize matrix.

Note that this method does not maintain data in the matrix.

Parameters
sizedimension

Definition at line 359 of file JMatrixND.hh.

360  {
361  static_cast<JMatrixND_t&>(*this).resize(size);
362 
363  getInstance().resize(size);
364  }
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:154
static data_type & getInstance()
Get unique instance of template class.
Definition: JSingleton.hh:27
void resize(const size_t size)
Resize matrix.
Definition: JMatrixND.hh:79
Basic NxN matrix.
Definition: JMatrixND.hh:35
JMatrixND& JMATH::JMatrixND::reset ( )
inlineinherited

Set matrix to the null matrix.

Returns
this matrix

Definition at line 372 of file JMatrixND.hh.

373  {
375 
376  double* p = A.data();
377 
378  for (size_t i = this->size(); i != 0; --i, ++p) {
379  *p = 0.0;
380  }
381 
382  p = this->data();
383 
384  for (size_t i = this->size(); i != 0; --i, p += this->size()) {
385  memcpy(p, A.data(), this->size() * sizeof(double));
386  }
387 
388  return *this;
389  }
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:154
static data_type & getInstance()
Get unique instance of template class.
Definition: JSingleton.hh:27
const double * data() const
Get pointer to data.
Definition: JMatrixND.hh:187
Basic NxN matrix.
Definition: JMatrixND.hh:35
source $JPP_DIR setenv csh $JPP_DIR eval JShellParser o a A
JMatrixND& JMATH::JMatrixND::setIdentity ( )
inlineinherited

Set to identity matrix.

Returns
this matrix

Definition at line 397 of file JMatrixND.hh.

398  {
399  reset();
400 
401  for (size_t i = 0; i != this->size(); ++i) {
402  (*this)(i,i) = 1.0;
403  }
404 
405  return *this;
406  }
JMatrixND & reset()
Set matrix to the null matrix.
Definition: JMatrixND.hh:372
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:154
JMatrixND& JMATH::JMatrixND::negate ( )
inlineinherited

Negate matrix.

Returns
this matrix

Definition at line 414 of file JMatrixND.hh.

415  {
416  double* p = this->data();
417 
418  for (size_t i = this->size()*this->size(); i != 0; --i, ++p) {
419  *p = -*p;
420  }
421 
422  return *this;
423  }
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:154
const double * data() const
Get pointer to data.
Definition: JMatrixND.hh:187
JMatrixND& JMATH::JMatrixND::add ( const JMatrixND A)
inlineinherited

Matrix addition.

Parameters
Amatrix
Returns
this matrix

Definition at line 432 of file JMatrixND.hh.

433  {
434  if (this->size() == A.size()) {
435 
436  double* p = this->data();
437  const double* q = A.data();
438 
439  for (size_t i = this->size()*this->size(); i != 0; --i, ++p, ++q) {
440  *p += *q;
441  }
442 
443  } else {
444  THROW(JException, "JMatrixND::add() inconsistent matrix dimensions " << this->size() << ' ' << A.size());
445  }
446  }
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:670
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:154
const double * data() const
Get pointer to data.
Definition: JMatrixND.hh:187
JMatrixND& JMATH::JMatrixND::sub ( const JMatrixND A)
inlineinherited

Matrix subtraction.

Parameters
Amatrix
Returns
this matrix

Definition at line 455 of file JMatrixND.hh.

456  {
457  if (this->size() == A.size()) {
458 
459  double* p = this->data();
460  const double* q = A.data();
461 
462  for (size_t i = this->size()*this->size(); i != 0; --i, ++p, ++q) {
463  *p -= *q;
464  }
465 
466  } else {
467  THROW(JException, "JMatrixND::sub() inconsistent matrix dimensions " << this->size() << ' ' << A.size());
468  }
469 
470  return *this;
471  }
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:670
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:154
const double * data() const
Get pointer to data.
Definition: JMatrixND.hh:187
JMatrixND& JMATH::JMatrixND::mul ( const double  factor)
inlineinherited

Scale matrix.

Parameters
factorfactor
Returns
this matrix

Definition at line 480 of file JMatrixND.hh.

481  {
482  double* p = this->data();
483 
484  for (size_t i = this->size()*this->size(); i != 0; --i, ++p) {
485  *p *= factor;
486  }
487 
488  return *this;
489  }
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:154
const double * data() const
Get pointer to data.
Definition: JMatrixND.hh:187
JMatrixND& JMATH::JMatrixND::mul ( const JMatrixND A,
const JMatrixND B 
)
inlineinherited

Matrix multiplication.

Parameters
Amatrix
Bmatrix
Returns
this matrix

Definition at line 517 of file JMatrixND.hh.

519  {
520  if (A.size() == B.size()) {
521 
522  this->resize(A.size());
523 
524  if (!this->empty()) {
525 
527 
528  C.set(B);
529  C.transpose();
530 
531  size_t i, row;
532 
533  for (row = 0; row + 4 <= this->size(); row += 4) { // process rows by four
534 
535  double* p0 = (*this)[row + 0];
536  double* p1 = (*this)[row + 1];
537  double* p2 = (*this)[row + 2];
538  double* p3 = (*this)[row + 3];
539 
540  for (size_t col = 0; col != this->size(); ++col, ++p0, ++p1, ++p2, ++p3) {
541 
542  double w0 = 0.0;
543  double w1 = 0.0;
544  double w2 = 0.0;
545  double w3 = 0.0;
546 
547  const double* a0 = A[row + 0];
548  const double* a1 = A[row + 1];
549  const double* a2 = A[row + 2];
550  const double* a3 = A[row + 3];
551 
552  const double* c0 = C[col];
553 
554  for (i = 0; i + 4 <= this->size(); i += 4, a0 += 4, a1 += 4, a2 += 4, a3 += 4, c0 += 4) {
555  w0 += a0[0] * c0[0] + a0[1] * c0[1] + a0[2] * c0[2] + a0[3] * c0[3];
556  w1 += a1[0] * c0[0] + a1[1] * c0[1] + a1[2] * c0[2] + a1[3] * c0[3];
557  w2 += a2[0] * c0[0] + a2[1] * c0[1] + a2[2] * c0[2] + a2[3] * c0[3];
558  w3 += a3[0] * c0[0] + a3[1] * c0[1] + a3[2] * c0[2] + a3[3] * c0[3];
559  }
560 
561  for ( ; i != this->size(); ++i, ++a0, ++a1, ++a2, ++a3, ++c0) {
562  w0 += (*a0) * (*c0);
563  w1 += (*a1) * (*c0);
564  w2 += (*a2) * (*c0);
565  w3 += (*a3) * (*c0);
566  }
567 
568  *p0 = w0;
569  *p1 = w1;
570  *p2 = w2;
571  *p3 = w3;
572  }
573  }
574 
575  for ( ; row != this->size(); ++row) { // remaining rows
576 
577  double* p0 = (*this)[row + 0];
578 
579  for (size_t col = 0; col != this->size(); ++col, ++p0) {
580 
581  double w0 = 0.0;
582 
583  const double* a0 = A[row + 0];
584  const double* c0 = C[col];
585 
586  for (i = 0; i + 4 <= this->size(); i += 4, a0 += 4, c0 += 4) {
587  w0 += a0[0] * c0[0] + a0[1] * c0[1] + a0[2] * c0[2] + a0[3] * c0[3];
588  }
589 
590  for ( ; i != this->size(); ++i, ++a0, ++c0) {
591  w0 += (*a0) * (*c0);
592  }
593 
594  *p0 = w0;
595  }
596  }
597  }
598 
599  } else {
600  THROW(JException, "JMatrixND::mul() inconsistent matrix dimensions " << A.size() << ' ' << B.size());
601  }
602 
603  return *this;
604  }
TPaveText * p1
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:670
void resize(const size_t size)
Resize matrix.
Definition: JMatrixND.hh:359
void set(const JMatrixND_t &A)
Set matrix.
Definition: JMatrixND.hh:107
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:154
static data_type & getInstance()
Get unique instance of template class.
Definition: JSingleton.hh:27
JMatrixND_t & transpose()
Transpose.
Definition: JMatrixND.hh:135
static const double C
Speed of light in vacuum [m/ns].
Definition: JConstants.hh:22
bool empty() const
Check emptyness.
Definition: JMatrixND.hh:176
Basic NxN matrix.
Definition: JMatrixND.hh:35
JMatrixND & JMATH::JMath< JMatrixND , JNullType >::mul ( const JNullType object)
inlineinherited

Multiply with object.

Parameters
objectobject
Returns
result object

Definition at line 273 of file JMath.hh.

274  {
275  return static_cast<JFirst_t&>(*this) = JCalculator<JFirst_t>::calculator.mul(static_cast<const JFirst_t&>(*this), object);
276  }
Auxiliary class for arithmetic operations on objects.
Definition: JCalculator.hh:18
JMatrixND& JMATH::JMatrixND::div ( const double  factor)
inlineinherited

Scale matrix.

Parameters
factorfactor
Returns
this matrix

Definition at line 498 of file JMatrixND.hh.

499  {
500  double* p = this->data();
501 
502  for (size_t i = this->size()*this->size(); i != 0; --i, ++p) {
503  *p /= factor;
504  }
505 
506  return *this;
507  }
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:154
const double * data() const
Get pointer to data.
Definition: JMatrixND.hh:187
bool JMATH::JMatrixND::equals ( const JMatrixND A,
const double  eps = std::numeric_limits<double>::min() 
) const
inlineinherited

Equality.

Parameters
Amatrix
epsnumerical precision
Returns
true if matrices identical; else false

Definition at line 614 of file JMatrixND.hh.

616  {
617  if (this->size() == A.size()) {
618 
619  for (size_t row = 0; row != this->size(); ++row) {
620 
621  const double* p = (*this)[row];
622  const double* q = A [row];
623 
624  for (size_t i = this->size(); i != 0; --i, ++p, ++q) {
625  if (fabs(*p - *q) > eps) {
626  return false;
627  }
628  }
629  }
630 
631  return true;
632 
633  } else {
634  THROW(JException, "JMatrixND::equals() inconsistent matrix dimensions " << this->size() << ' ' << A.size());
635  }
636  }
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:670
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:154
bool JMATH::JMatrixND::isIdentity ( const double  eps = std::numeric_limits<double>::min()) const
inlineinherited

Test identity.

Parameters
epsnumerical precision
Returns
true if identity matrix; else false

Definition at line 645 of file JMatrixND.hh.

646  {
647  for (size_t i = 0; i != this->size(); ++i) {
648 
649  if (fabs(1.0 - (*this)(i,i)) > eps) {
650  return false;
651  };
652 
653  for (size_t j = 0; j != i; ++j) {
654  if (fabs((*this)(i,j)) > eps || fabs((*this)(j,i)) > eps) {
655  return false;
656  };
657  }
658  }
659 
660  return true;
661  }
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:154
int j
Definition: JPolint.hh:634
bool JMATH::JMatrixND::isSymmetric ( const double  eps = std::numeric_limits<double>::min()) const
inlineinherited

Test symmetry.

Parameters
epsnumerical precision
Returns
true if symmetric matrix; else false

Definition at line 670 of file JMatrixND.hh.

671  {
672  for (size_t i = 0; i != this->size(); ++i) {
673  for (size_t j = 0; j != i; ++j) {
674  if (fabs((*this)(i,j) - (*this)(j,i)) > eps) {
675  return false;
676  };
677  }
678  }
679 
680  return true;
681  }
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:154
int j
Definition: JPolint.hh:634
double JMATH::JMatrixND::getDot ( const JVectorND v) const
inlineinherited

Get dot product.

The dot product corresponds to

\[ v^{T} \times A \times v \]

.

Parameters
vvector
Returns
dot product

Definition at line 692 of file JMatrixND.hh.

693  {
694  double dot = 0.0;
695 
696  for (size_t i = 0; i != v.size(); ++i) {
697 
698  const double* p = (*this)[i];
699 
700  double w = 0.0;
701 
702  for (JVectorND::const_iterator y = v.begin(); y != v.end(); ++y, ++p) {
703  w += (*p) * (*y);
704  }
705 
706  dot += v[i] * w;
707  }
708 
709  return dot;
710  }
data_type w[N+1][M+1]
Definition: JPolint.hh:708
void JMATH::JMatrixND::rswap ( size_t  r1,
size_t  r2 
)
inlineprotectedinherited

Definition at line 792 of file JMatrixND.hh.

793  {
795 
796  memcpy(A.data(), (*this)[r1], this->size() * sizeof(double));
797  memcpy((*this)[r1], (*this)[r2], this->size() * sizeof(double));
798  memcpy((*this)[r2], A.data(), this->size() * sizeof(double));
799  }
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:154
static data_type & getInstance()
Get unique instance of template class.
Definition: JSingleton.hh:27
const double * data() const
Get pointer to data.
Definition: JMatrixND.hh:187
Basic NxN matrix.
Definition: JMatrixND.hh:35
source $JPP_DIR setenv csh $JPP_DIR eval JShellParser o a A
void JMATH::JMatrixND::cswap ( size_t  c1,
size_t  c2 
)
inlineprotectedinherited

Swap columns.

Parameters
c1column
c2column

Definition at line 808 of file JMatrixND.hh.

809  {
810  using std::swap;
811 
812  double* p1 = this->data() + c1;
813  double* p2 = this->data() + c2;
814 
815  for (size_t i = this->size(); i != 0; --i, p1 += this->size(), p2 += this->size()) {
816  swap(*p1, *p2);
817  }
818  }
TPaveText * p1
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:154
TCanvas * c1
Global variables to handle mouse events.
void swap(JMatrixND_t &A)
Swap matrices.
Definition: JMatrixND.hh:120
const double * data() const
Get pointer to data.
Definition: JMatrixND.hh:187
void JMATH::JMatrixND_t::clear ( )
inlineinherited

Clear memory.

Definition at line 60 of file JMatrixND.hh.

61  {
62  if (__p != NULL) {
63  delete [] __p;
64  }
65 
66  __p = NULL;
67  __n = 0;
68  __m = 0;
69  }
size_t __n
dimension of matrix
Definition: JMatrixND.hh:292
double * __p
pointer to data
Definition: JMatrixND.hh:291
size_t __m
capacity of matrix
Definition: JMatrixND.hh:293
void JMATH::JMatrixND_t::set ( const JMatrixND_t A)
inlineinherited

Set matrix.

Parameters
Amatrix

Definition at line 107 of file JMatrixND.hh.

108  {
109  this->resize(A.size());
110 
111  memcpy(this->data(), A.data(), A.size() * A.size() * sizeof(double));
112  }
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:154
void resize(const size_t size)
Resize matrix.
Definition: JMatrixND.hh:79
const double * data() const
Get pointer to data.
Definition: JMatrixND.hh:187
void JMATH::JMatrixND_t::swap ( JMatrixND_t A)
inlineinherited

Swap matrices.

Parameters
Amatrix

Definition at line 120 of file JMatrixND.hh.

121  {
122  using std::swap;
123 
124  swap(this->__p, A.__p);
125  swap(this->__n, A.__n);
126  swap(this->__m, A.__m);
127  }
size_t __n
dimension of matrix
Definition: JMatrixND.hh:292
double * __p
pointer to data
Definition: JMatrixND.hh:291
void swap(JMatrixND_t &A)
Swap matrices.
Definition: JMatrixND.hh:120
size_t __m
capacity of matrix
Definition: JMatrixND.hh:293
JMatrixND_t& JMATH::JMatrixND_t::transpose ( )
inlineinherited

Transpose.

Returns
this matrix

Definition at line 135 of file JMatrixND.hh.

136  {
137  using std::swap;
138 
139  for (size_t row = 0; row != this->size(); ++row) {
140  for (size_t col = 0; col != row; ++col) {
141  swap((*this)(row,col), (*this)(col,row));
142  }
143  }
144 
145  return *this;
146  }
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:154
void swap(JMatrixND_t &A)
Swap matrices.
Definition: JMatrixND.hh:120
size_t JMATH::JMatrixND_t::size ( ) const
inlineinherited

Get dimension of matrix.

Returns
dimension

Definition at line 154 of file JMatrixND.hh.

155  {
156  return __n;
157  }
size_t __n
dimension of matrix
Definition: JMatrixND.hh:292
size_t JMATH::JMatrixND_t::capacity ( ) const
inlineinherited

Get capacity of dimension.

Returns
capacity

Definition at line 165 of file JMatrixND.hh.

166  {
167  return __m;
168  }
size_t __m
capacity of matrix
Definition: JMatrixND.hh:293
bool JMATH::JMatrixND_t::empty ( ) const
inlineinherited

Check emptyness.

Returns
true if empty; else false

Definition at line 176 of file JMatrixND.hh.

177  {
178  return __n == 0;
179  }
size_t __n
dimension of matrix
Definition: JMatrixND.hh:292
const double* JMATH::JMatrixND_t::data ( ) const
inlineinherited

Get pointer to data.

Returns
pointer to data.

Definition at line 187 of file JMatrixND.hh.

188  {
189  return __p;
190  }
double * __p
pointer to data
Definition: JMatrixND.hh:291
double* JMATH::JMatrixND_t::data ( )
inlineinherited

Get pointer to data.

Returns
pointer to data.

Definition at line 198 of file JMatrixND.hh.

199  {
200  return __p;
201  }
double * __p
pointer to data
Definition: JMatrixND.hh:291
const double* JMATH::JMatrixND_t::operator[] ( size_t  row) const
inlineinherited

Get row data.

Parameters
rowrow number
Returns
pointer to row data

Definition at line 210 of file JMatrixND.hh.

211  {
212  return __p + row*__n;
213  }
size_t __n
dimension of matrix
Definition: JMatrixND.hh:292
double * __p
pointer to data
Definition: JMatrixND.hh:291
double* JMATH::JMatrixND_t::operator[] ( size_t  row)
inlineinherited

Get row data.

Parameters
rowrow number
Returns
pointer to row data

Definition at line 222 of file JMatrixND.hh.

223  {
224  return __p + row*__n;
225  }
size_t __n
dimension of matrix
Definition: JMatrixND.hh:292
double * __p
pointer to data
Definition: JMatrixND.hh:291
double JMATH::JMatrixND_t::operator() ( const size_t  row,
const size_t  col 
) const
inlineinherited

Get matrix element.

Parameters
rowrow number
colcolumn number
Returns
matrix element at (row,col)

Definition at line 235 of file JMatrixND.hh.

236  {
237  return *(__p + row*__n + col);
238  }
size_t __n
dimension of matrix
Definition: JMatrixND.hh:292
double * __p
pointer to data
Definition: JMatrixND.hh:291
double& JMATH::JMatrixND_t::operator() ( const size_t  row,
const size_t  col 
)
inlineinherited

Get matrix element.

Parameters
rowrow number
colcolumn number
Returns
matrix element at (row,col)

Definition at line 248 of file JMatrixND.hh.

249  {
250  return *(__p + row*__n + col);
251  }
size_t __n
dimension of matrix
Definition: JMatrixND.hh:292
double * __p
pointer to data
Definition: JMatrixND.hh:291
double JMATH::JMatrixND_t::at ( size_t  row,
size_t  col 
) const
inlineinherited

Get matrix element.

Parameters
rowrow number
colcolumn number
Returns
matrix element at (row,col)

Definition at line 261 of file JMatrixND.hh.

262  {
263  if (row >= this->size() ||
264  col >= this->size()) {
265  THROW(JIndexOutOfRange, "JMatrixND::at(" << row << "," << col << "): index out of range.");
266  }
267 
268  return (*this)(row, col);
269  }
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:670
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:154
double& JMATH::JMatrixND_t::at ( size_t  row,
size_t  col 
)
inlineinherited

Get matrix element.

Parameters
rowrow number
colcolumn number
Returns
matrix element at (row,col)

Definition at line 279 of file JMatrixND.hh.

280  {
281  if (row >= this->size() ||
282  col >= this->size()) {
283  THROW(JIndexOutOfRange, "JMatrixND::at(" << row << "," << col << "): index out of range.");
284  }
285 
286  return (*this)(row, col);
287  }
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:670
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:154
template<class T>
static data_type& JLANG::JSingleton< T >::getInstance ( )
inlinestaticinherited

Get unique instance of template class.

Returns
object

Definition at line 27 of file JSingleton.hh.

28  {
29  static data_type value;
30 
31  return value;
32  }

Member Data Documentation

std::vector<variance> JFIT::JMatrixNZ::buffer
private

Definition at line 186 of file JMatrixNZ.hh.

double* JMATH::JMatrixND_t::__p
protectedinherited

pointer to data

Definition at line 291 of file JMatrixND.hh.

size_t JMATH::JMatrixND_t::__n
protectedinherited

dimension of matrix

Definition at line 292 of file JMatrixND.hh.

size_t JMATH::JMatrixND_t::__m
protectedinherited

capacity of matrix

Definition at line 293 of file JMatrixND.hh.


The documentation for this class was generated from the following file: