Belle II Software  release-06-02-00
kinem.h
1 /*
2  * kinem.h - kinematics classes
3  *
4  * this file is part of PJFry library
5  * Copyright 2011 Valery Yundin
6  */
7 
8 #ifndef QUL_KINEM_H
9 #define QUL_KINEM_H
10 
11 #include "common.h"
12 #include <cstdlib>
13 // #include <cfloat>
14 
15 template <int R>
16 class Kinem {
17 public:
18  bool operator == (const Kinem& kinem) const;
19 
20  inline
21  double mass(int i) const
22  {
23  return kdata[(i - 1) * (i + 2) / 2];
24  }
25 
26 protected:
27  Kinem() {}
28 
29  void zero(); // initialize to zero
30 
31  typedef union { double d64; int64_t i64; } DI64;
32 // static const double mdelta=1e-14;
33 // static const int64_t idelta=1+mdelta/DBL_EPSILON; // needs cfloat
34  /*
35  Comparison tolerance for double
36  idelta fdelta(%)
37  5 - 1.11e-15
38  46 - 1.02e-14
39  451 - 1.00e-13
40  4504 - 1.00e-12
41  45036 - 1.00e-11
42  */
43  static const uint64_t idelta = 5; // unsigned !!!
44 
45  inline
46  bool eq(const double& a, const double& b) const
47  {
48  const DI64 ia = {a};
49  const DI64 ib = {b};
50  const int64_t diff = ia.i64 - ib.i64;
51  return diff == 0LL || static_cast<uint64_t>(llabs(diff)) <= idelta;
52  }
53  inline
54  bool neq(const double& a, const double& b) const
55  {
56  const DI64 ia = {a};
57  const DI64 ib = {b};
58  const int64_t diff = ia.i64 - ib.i64;
59  return diff != 0LL && static_cast<uint64_t>(llabs(diff)) > idelta;
60  }
61 
62  static const int KLEN = R * (R + 1) / 2;
63  enum Invar {im1 = 0,
64  ip1, im2,
65  ip2, ip3, im3,
66  ip4, is12, is23, im4,
67  ip5, is34, is45, is15, im5,
68  ip6, /* 6-point inv */ im6
69  };
70 
71  double kdata[KLEN];
72 };
73 
74 template <int R>
75 bool Kinem<R>::operator == (const Kinem<R>& kinem) const
76 {
77  for (int i = 0; i < KLEN; i++) {
78  if (not eq(kdata[i], kinem.kdata[i])) return false;
79  }
80  return true;
81 }
82 
83 template <int R>
84 void Kinem<R>::zero()
85 {
86  for (int i = 0; i < KLEN; i++) {
87  kdata[i] = 0;
88  }
89 }
90 
91 // 1-point kinematics
92 class Kinem1 : public Kinem<1> {
93 public:
94  Kinem1() { zero(); }
95  Kinem1(double xm1)
96  {
97  kdata[im1] = xm1;
98  }
99 
100  inline double m1() const { return kdata[im1]; }
101 };
102 
103 // 2-point kinematics
104 class Kinem2 : public Kinem<2> {
105 public:
106  Kinem2() { zero(); }
107  Kinem2(double xp1,
108  double xm1, double xm2)
109  {
110  kdata[ip1] = xp1;
111  kdata[im1] = xm1;
112  kdata[im2] = xm2;
113  }
114 
115  inline double p1() const { return kdata[ip1]; }
116  inline double m1() const { return kdata[im1]; }
117  inline double m2() const { return kdata[im2]; }
118 };
119 
120 // 3-point kinematics
121 class Kinem3 : public Kinem<3> {
122 public:
123  Kinem3() { zero(); }
124  Kinem3(double xp1, double xp2, double xp3,
125  double xm1, double xm2, double xm3)
126  {
127  kdata[ip1] = xp1;
128  kdata[ip2] = xp2;
129  kdata[ip3] = xp3;
130  kdata[im1] = xm1;
131  kdata[im2] = xm2;
132  kdata[im3] = xm3;
133  }
134 
135  inline double p1() const { return kdata[ip1]; }
136  inline double p2() const { return kdata[ip2]; }
137  inline double p3() const { return kdata[ip3]; }
138  inline double m1() const { return kdata[im1]; }
139  inline double m2() const { return kdata[im2]; }
140  inline double m3() const { return kdata[im3]; }
141 };
142 
143 // 4-point kinematics
144 class Kinem4 : public Kinem<4> {
145 public:
146  Kinem4() { zero(); }
147  Kinem4(double xp1, double xp2, double xp3, double xp4,
148  double xs12, double xs23,
149  double xm1, double xm2, double xm3, double xm4)
150  {
151  kdata[ip1] = xp1;
152  kdata[ip2] = xp2;
153  kdata[ip3] = xp3;
154  kdata[ip4] = xp4;
155  kdata[is12] = xs12;
156  kdata[is23] = xs23;
157  kdata[im1] = xm1;
158  kdata[im2] = xm2;
159  kdata[im3] = xm3;
160  kdata[im4] = xm4;
161  }
162 
163  inline double p1() const { return kdata[ip1]; }
164  inline double p2() const { return kdata[ip2]; }
165  inline double p3() const { return kdata[ip3]; }
166  inline double p4() const { return kdata[ip4]; }
167  inline double m1() const { return kdata[im1]; }
168  inline double m2() const { return kdata[im2]; }
169  inline double m3() const { return kdata[im3]; }
170  inline double m4() const { return kdata[im4]; }
171  inline double s12() const { return kdata[is12]; }
172  inline double s23() const { return kdata[is23]; }
173 };
174 
175 // 5-point kinematics
176 class Kinem5 : public Kinem<5> {
177 public:
178  Kinem5() { zero(); }
179  Kinem5(double xp1, double xp2, double xp3, double xp4, double xp5,
180  double xs12, double xs23, double xs34, double xs45, double xs15,
181  double xm1, double xm2, double xm3, double xm4, double xm5)
182  {
183  kdata[ip1] = xp1;
184  kdata[ip2] = xp2;
185  kdata[ip3] = xp3;
186  kdata[ip4] = xp4;
187  kdata[ip5] = xp5;
188  kdata[is12] = xs12;
189  kdata[is23] = xs23;
190  kdata[is34] = xs34;
191  kdata[is45] = xs45;
192  kdata[is15] = xs15;
193  kdata[im1] = xm1;
194  kdata[im2] = xm2;
195  kdata[im3] = xm3;
196  kdata[im4] = xm4;
197  kdata[im5] = xm5;
198  }
199 
200  inline double p1() const { return kdata[ip1]; }
201  inline double p2() const { return kdata[ip2]; }
202  inline double p3() const { return kdata[ip3]; }
203  inline double p4() const { return kdata[ip4]; }
204  inline double p5() const { return kdata[ip5]; }
205  inline double m1() const { return kdata[im1]; }
206  inline double m2() const { return kdata[im2]; }
207  inline double m3() const { return kdata[im3]; }
208  inline double m4() const { return kdata[im4]; }
209  inline double m5() const { return kdata[im5]; }
210  inline double s12() const { return kdata[is12]; }
211  inline double s23() const { return kdata[is23]; }
212  inline double s34() const { return kdata[is34]; }
213  inline double s45() const { return kdata[is45]; }
214  inline double s15() const { return kdata[is15]; }
215 };
216 
217 #endif /* QUL_KINEM_H */
Kinem1
Definition: kinem.h:92
Kinem2
Definition: kinem.h:104
Kinem3
Definition: kinem.h:121
Kinem4
Definition: kinem.h:144
Kinem5
Definition: kinem.h:176
Kinem
Definition: kinem.h:16