Jpp 19.3.0-rc.2
the software that should make you happy
Loading...
Searching...
No Matches
JCompass/JModel.hh
Go to the documentation of this file.
1#ifndef __JCOMPASS_JMODEL__
2#define __JCOMPASS_JMODEL__
3
6#include "JLang/JException.hh"
7#include "JMath/JMath.hh"
10#include "JFit/JMEstimator.hh"
11
12#include "JCompass/JHit.hh"
13
14
15/**
16 * \author mdejong
17 */
18namespace JCOMPASS {}
19namespace JPP { using namespace JCOMPASS; }
20
21namespace JCOMPASS {
22
23 using JMATH::JMath;
28
29 /**
30 * Model.
31 */
32 struct JModel :
33 public JMath<JModel>
34 {
35 /**
36 * Default constructor.
37 */
39 Q0(JQuaternion3D::getIdentity()),
40 Q1(JQuaternion3D::getIdentity())
41 {}
42
43
44 /**
45 * Constructor.
46 *
47 * \param Q0 tilt
48 * \param Q1 twist
49 */
51 const JQuaternion3D& Q1) :
52 Q0(Q0),
53 Q1(Q1)
54 {}
55
56
57 /**
58 * Constructor.
59 *
60 * The data type corresponding to the hits should provide for the following policy methods.
61 * <pre>
62 * double getZ(); // get z-position
63 * JQuaternion3D getQuaternion(); // get quaternion
64 * </pre>
65 * Note that the input data should be ordered with increasing z-positions.
66 *
67 * \param __begin begin of data
68 * \param __end end of data
69 */
70 template<class T>
71 JModel(T __begin,
72 T __end) :
73 Q0(JQuaternion3D::getIdentity()),
74 Q1(JQuaternion3D::getIdentity())
75 {
76 using namespace std;
77 using namespace JPP;
78
79 const int N = distance(__begin, __end);
80
81 if (N >= NUMBER_OF_PARAMETERS) {
82
84
85 for (T q = __begin, p = q++; q != __end; ++p, ++q) {
86
87 const double dz = q->getZ() - p->getZ();
88
89 if (dz != 0.0) {
90
91 JQuaternion3D Q(p->getQuaternion());
92
93 Q.conjugate();
94 Q.mul(q->getQuaternion());
95 Q.pow(1.0 / dz);
96
97 buffer.push_back(Q);
98 }
99 }
100
101 Q1 = getAverage(buffer.begin(), buffer.end());
102 Q1 = JQuaternion3D::decomposition(Q1, JVector3Z_t).twist;
103
104 const double z1 = getAverage(make_array(__begin, __end, &JHit::getZ));
105
106 Q0 = getAverage(make_array(__begin, __end, &JHit::getQuaternion));
107 Q0 = pow(Q1, -z1) * Q0;
108
109 } else {
110 throw JValueOutOfRange("JModel: Not enough data points.");
111 }
112 }
113
114
115 /**
116 * Add model.
117 *
118 * \param model model
119 * \return this model
120 */
121 JModel& add(const JModel& model)
122 {
123 Q0 *= model.Q0;
124 Q1 *= model.Q1;
125
126 Q0.normalise();
127 Q1.normalise();
128
129 return *this;
130 }
131
132
133 /**
134 * Subtract model.
135 *
136 * \param model model
137 * \return this model
138 */
139 JModel& sub(const JModel& model)
140 {
141 Q0 *= model.Q0.getConjugate();
142 Q1 *= model.Q1.getConjugate();
143
144 Q0.normalise();
145 Q1.normalise();
146
147 return *this;
148 }
149
150
151 /**
152 * Scale model.
153 *
154 * \param factor multiplication factor
155 * \return this model
156 */
157 JModel& mul(const double factor)
158 {
159 Q0.pow(factor);
160 Q1.pow(factor);
161
162 return *this;
163 }
164
165
166 /**
167 * Get quaternion at given z-position.
168 *
169 * \param z z-position.
170 * \return quaternion
171 */
172 JQuaternion3D operator()(const double z) const
173 {
174 using namespace JPP;
175
176 return pow(Q1, z) * Q0;
177 }
178
179 static const int NUMBER_OF_PARAMETERS = 4; //!< number of parameters of fit per quaternion
180
181 JQuaternion3D Q0; //!< tilt
182 JQuaternion3D Q1; //!< twist
183 };
184
185
186 /**
187 * Auxiliary data structure for chi2 evaluation.
188 */
189 struct JChi2 {
190 /**
191 * Constructor.
192 *
193 * \param type M-Estimator type
194 */
195 JChi2(const int type) :
196 estimator(getMEstimator(type))
197 {}
198
199
200 /**
201 * Fit function.
202 *
203 * \param model model
204 * \param hit hit
205 * \return chi2
206 */
207 inline double operator()(const JModel& model, const JHit& hit) const
208 {
209 return estimator->getRho(getAngle(model(hit.getZ()), hit.getQuaternion()) / hit.getSigma());
210 }
211
213 };
214}
215
216#endif
Exceptions.
Maximum likelihood estimator (M-estimators).
Auxiliary methods for geometrical methods.
Base class for data structures with artithmetic capabilities.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Data structure for unit quaternion in three dimensions.
const JQuaternion3D & getQuaternion() const
Get quaternion.
JQuaternion3D & normalise()
Normalise quaternion.
JQuaternion3D & conjugate()
Conjugate quaternion.
JQuaternion3D & mul(const double factor)
Scale quaternion.
JQuaternion3D & pow(const double y)
Raise quaternion to given power.
The template JSharedPointer class can be used to share a pointer to an object.
Exception for accessing a value in a collection that is outside of its range.
Auxiliary classes and methods for orientation calibration based on compasses.
JMEstimator * getMEstimator(const int type)
Get M-Estimator.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Auxiliary data structure for chi2 evaluation.
double operator()(const JModel &model, const JHit &hit) const
Fit function.
JChi2(const int type)
Constructor.
JLANG::JSharedPointer< JMEstimator > estimator
M-Estimator function.
double getSigma() const
Get resolution.
double getZ() const
Get z-position.
JModel & sub(const JModel &model)
Subtract model.
JModel & add(const JModel &model)
Add model.
JQuaternion3D operator()(const double z) const
Get quaternion at given z-position.
JQuaternion3D Q0
tilt
JModel(T __begin, T __end)
Constructor.
static const int NUMBER_OF_PARAMETERS
number of parameters of fit per quaternion
JQuaternion3D Q1
twist
JModel & mul(const double factor)
Scale model.
JModel()
Default constructor.
JModel(const JQuaternion3D &Q0, const JQuaternion3D &Q1)
Constructor.
Interface for maximum likelihood estimator (M-estimator).
Auxiliary data structure for decomposition of quaternion in twist and swing quaternions.
JQuaternion3D twist
rotation around parallel axis
Auxiliary base class for aritmetic operations of derived class types.
Definition JMath.hh:347