Jpp  18.0.0-rc.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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"
8 #include "JMath/JMathToolkit.hh"
10 #include "JFit/JMEstimator.hh"
11 
12 #include "JCompass/JHit.hh"
13 
14 
15 /**
16  * \author mdejong
17  */
18 namespace JCOMPASS {}
19 namespace JPP { using namespace JCOMPASS; }
20 
21 namespace JCOMPASS {
22 
23  using JMATH::JMath;
26  using JFIT::JMEstimator;
27  using JFIT::getMEstimator;
28 
29  /**
30  * Model.
31  */
32  struct JModel :
33  public JMath<JModel>
34  {
35  /**
36  * Default constructor.
37  */
38  JModel() :
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 
83  vector<JQuaternion3D> buffer;
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());
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
double getAngle(const JQuaternion3D &first, const JQuaternion3D &second)
Get space angle between quanternions.
JQuaternion3D & pow(const double y)
Raise quaternion to given power.
Exceptions.
Q(UTCMax_s-UTCMin_s)-livetime_s
Auxiliary base class for aritmetic operations of derived class types.
Definition: JMath.hh:109
Auxiliary methods for geometrical methods.
JQuaternion3D getConjugate() const
Get conjugate of this quaternion.
Interface for maximum likelihood estimator (M-estimator).
Definition: JMEstimator.hh:20
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
JLANG::JSharedPointer< JMEstimator > estimator
M-Estimator function.
JModel & add(const JModel &model)
Add model.
std::iterator_traits< T >::value_type getAverage(T __begin, T __end)
Get average.
Definition: JMath.hh:494
double getZ() const
Get z-position.
const JQuaternion3D & getQuaternion() const
Get quaternion.
Auxiliary data structure for chi2 evaluation.
The template JSharedPointer class can be used to share a pointer to an object.
JQuaternion3D twist
rotation around parallel axis
JQuaternion3D Q0
tilt
const array_type< JValue_t > & make_array(const JValue_t(&array)[N])
Method to create array of values.
Definition: JVectorize.hh:54
do set_variable OUTPUT_DIRECTORY $WORKDIR T
T pow(const T &x, const double y)
Power .
Definition: JMath.hh:97
JModel()
Default constructor.
JModel & mul(const double factor)
Scale model.
JChi2(const int type)
Constructor.
JModel(const JQuaternion3D &Q0, const JQuaternion3D &Q1)
Constructor.
Data structure for unit quaternion in three dimensions.
JQuaternion3D Q1
twist
then usage $script< input file >[option[primary[working directory]]] nWhere option can be N
Definition: JMuonPostfit.sh:36
static const JVector3D JVector3Z_t(0, 0, 1)
unit z-vector
Base class for data structures with artithmetic capabilities.
double operator()(const JModel &model, const JHit &hit) const
Fit function.
JQuaternion3D & normalise()
Normalise quaternion.
Exception for accessing a value in a collection that is outside of its range.
Definition: JException.hh:162
JModel(T __begin, T __end)
Constructor.
JMEstimator * getMEstimator(const int type)
Get M-Estimator.
Definition: JMEstimator.hh:203
JQuaternion3D operator()(const double z) const
Get quaternion at given z-position.
double getSigma() const
Get resolution.
static const int NUMBER_OF_PARAMETERS
number of parameters of fit per quaternion
Auxiliary data structure for decomposition of quaternion in twist and swing quaternions.
Maximum likelihood estimator (M-estimators).
JModel & sub(const JModel &model)
Subtract model.