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