Jpp  19.1.0-rc.1
the software that should make you happy
JMuonPrefit.hh
Go to the documentation of this file.
1 #ifndef __JRECONSTRUCTION__JMUONPREFIT__
2 #define __JRECONSTRUCTION__JMUONPREFIT__
3 
4 #include <iostream>
5 #include <iomanip>
6 #include <vector>
7 #include <algorithm>
8 
12 
13 #include "JTrigger/JHitR1.hh"
14 #include "JTrigger/JBuildL0.hh"
15 #include "JTrigger/JBuildL2.hh"
16 #include "JTrigger/JAlgorithm.hh"
17 #include "JTrigger/JMatch1D.hh"
18 #include "JTrigger/JMatch3B.hh"
19 
20 #include "JFit/JLine1Z.hh"
21 #include "JFit/JLine1ZEstimator.hh"
22 #include "JFit/JMatrixNZ.hh"
23 #include "JFit/JVectorNZ.hh"
24 #include "JFit/JFitToolkit.hh"
25 
26 #include "JReconstruction/JEvt.hh"
29 
30 #include "JMath/JConstants.hh"
31 #include "JTools/JPermutation.hh"
32 
33 #include "JLang/JPredicate.hh"
34 
37 
38 #include "JGeometry3D/JOmega3D.hh"
40 
41 #include "Jeep/JMessage.hh"
42 
43 
44 /**
45  * \author mdejong, gmaggi, azegarelli
46  */
47 
48 namespace JRECONSTRUCTION {}
49 namespace JPP { using namespace JRECONSTRUCTION; }
50 
51 namespace JRECONSTRUCTION {
52 
54  using JFIT::JLine1Z;
55  using JFIT::JEstimator;
56  using JFIT::JMatrixNZ;
57  using JFIT::JVectorNZ;
58  using JTRIGGER::JMatch3B;
59  using JTRIGGER::JMatch1D;
61 
62 
63  /**
64  * Wrapper class to make pre-fit of muon trajectory.
65  *
66  * The JMuonPrefit fit is used to generate start values for subsequent fits (usually JMuonSimplex and JMuonGandalf).\n
67  * To this end, a scan of directions is made and the time and transverse positions of the track are fitted for each direction (JFIT::JEstimator<JLine1Z>).\n
68  * The directions are spaced by the parameters JMuonPrefitParameters_t::gridAngle_deg.\n
69  * This angle corresponds to the envisaged angular accuracy of the result.\n
70  * The probability that one of the results is less than this angle away from the correct value,
71  * multiple start values should be considered (JMuonPrefitParameters_t::numberOfPrefits).\n
72  * Note that the CPU time scales with the inverse of the square of this angle.\n
73  * The chi-squared is based on the time residuals.\n
74  */
75  struct JMuonPrefit :
77  public JEstimator<JLine1Z>
78  {
82 
83  using JEstimator_t::operator();
84 
85  /**
86  * Constructor
87  *
88  * \param parameters parameters
89  * \param router module router
90  * \param debug debug
91  */
93  const JModuleRouter& router,
94  const int debug = 0) :
95  JMuonPrefitParameters_t(parameters),
96  router(router),
97  omega (parameters.gridAngle_deg * JMATH::PI/180.0),
98  debug (debug)
99  {
100  configure();
101  }
102 
103 
104  /**
105  * Constructor
106  *
107  * \param router module router
108  * \param parameters parameters
109  * \param omega directions
110  * \param debug debug
111  */
113  const JModuleRouter& router,
114  const JOmega3D& omega,
115  const int debug = 0) :
116  JMuonPrefitParameters_t(parameters),
117  router(router),
118  omega (omega),
119  debug (debug)
120  {
121  configure();
122  }
123 
124 
125  /**
126  * Fit function.
127  *
128  * \param event event
129  * \return fit results
130  */
132  {
133  using namespace std;
134  using namespace JPP;
135 
136  const JBuildL0<hit_type> buildL0;
138 
139  buffer_type dataL0;
140  buffer_type dataL1;
141 
142  buildL2(event, router, !useL0, back_inserter(dataL1));
143 
144  // 3D cluster of unique optical modules
145 
147 
148  sort(dataL1.begin(), dataL1.end(), hit_type::compare);
149 
150  buffer_type::iterator __end = dataL1.end();
151 
152  __end = unique(dataL1.begin(), __end, equal_to<JDAQModuleIdentifier>());
153 
154  __end = clusterizeWeight(dataL1.begin(), __end, match3B);
155 
156  dataL1.erase(__end, dataL1.end());
157 
158 
159  if (useL0) {
160 
161  buildL0(event, router, true, back_inserter(dataL0));
162 
163  __end = dataL0.end();
164 
165  for (buffer_type::iterator i = dataL0.begin(); i != __end; ) {
166 
167  if (match3B.count(*i, dataL1.begin(), dataL1.end()) != 0)
168  ++i;
169  else
170  swap(*i, *--__end);
171  }
172 
173  dataL0.erase(__end, dataL0.end());
174  }
175 
176  return (*this)(dataL0, dataL1);
177  }
178 
179 
180  /**
181  * Fit function.
182  *
183  * \param dataL0 L0 hit data
184  * \param dataL1 L1 hit data
185  * \return fit results
186  */
187  JEvt operator()(const buffer_type& dataL0,
188  const buffer_type& dataL1)
189  {
190  using namespace std;
191  using namespace JPP;
192 
193  const double STANDARD_DEVIATIONS = 3.0; // [unit]
194  const double HIT_OFF = 1.0e3 * sigma_ns * sigma_ns; // [ns^2]
195 
197 
198  data.reserve(dataL0.size() +
199  dataL1.size());
200 
201  JEvt out;
202 
203  for (JOmega3D_t::const_iterator dir = omega.begin(); dir != omega.end(); ++dir) {
204 
205  const JRotation3D R(*dir);
206 
207 
208  buffer_type::iterator __end = copy(dataL1.begin(), dataL1.end(), data.begin());
209 
210  for (buffer_type::iterator i = data.begin(); i != __end; ++i) {
211  i->rotate(R);
212  }
213 
214 
215  // reduce data
216 
217  if (distance(data.begin(), __end) > NMaxHits) {
218 
219  buffer_type::iterator __p = data.begin();
220 
221  advance(__p, NMaxHits);
222 
223  partial_sort(data.begin(), __p, __end, cmz);
224 
225  __end = __p;
226  }
227 
228 
229  // 1D cluster
230 
231  __end = clusterizeWeight(data.begin(), __end, match1D);
232 
233  if (useL0) {
234 
235  buffer_type::iterator p = __end; // begin L0 data
236  buffer_type::iterator q = copy(dataL0.begin(), dataL0.end(), p); // end L0 data
237 
238  for (buffer_type::iterator i = p; i != q; ++i) {
239 
240  if (find_if(data.begin(), __end, make_predicate(&hit_type::getModuleID, i->getModuleID())) == __end) {
241 
242  i->rotate(R);
243 
244  if (match1D.count(*i, data.begin(), __end) != 0) {
245  *p = *i;
246  ++p;
247  }
248  }
249  }
250 
251  __end = clusterize(__end, p, match1D);
252  }
253 
254 
255  if (distance(data.begin(), __end) <= NUMBER_OF_PARAMETERS) {
256  continue;
257  }
258 
259 
260  // 1D fit
261 
262  JLine1Z tz;
263  double chi2 = numeric_limits<double>::max();
264  int NDF = distance(data.begin(), __end) - NUMBER_OF_PARAMETERS;
265  int N = getCount(data.begin(), __end);
266 
267 
268  if (distance(data.begin(), __end) <= factoryLimit) {
269 
270  int number_of_outliers = numberOfOutliers;
271 
272  if (number_of_outliers > NDF - 1) {
273  number_of_outliers = NDF - 1;
274  }
275 
276  double ymin = numeric_limits<double>::max();
277 
278  buffer_type::iterator __end1 = __end;
279 
280  for (int n = 0; n <= number_of_outliers; ++n, --__end1) {
281 
282  sort(data.begin(), __end, hit_type::compare);
283 
284  do {
285  /*
286  if (getNumberOfStrings(router, data.begin(), __end1) < 2) {
287  continue;
288  }
289  */
290  try {
291 
292  (*this)(data.begin(), __end1);
293 
294  V.set(*this, data.begin(), __end1, gridAngle_deg, sigma_ns);
295  Y.set(*this, data.begin(), __end1);
296 
297  V.invert();
298 
299  double y = getChi2(Y, V);
300 
301  if (y <= -(STANDARD_DEVIATIONS * STANDARD_DEVIATIONS)) {
302 
303  WARNING(endl << "chi2(1) " << y << endl);
304 
305  } else {
306 
307  if (y < 0.0) {
308  y = 0.0;
309  }
310 
311  if (y < ymin) {
312  ymin = y;
313  tz = *this;
314  chi2 = ymin;
315  NDF = distance(data.begin(), __end1) - NUMBER_OF_PARAMETERS;
316  N = getCount(data.begin(), __end1);
317  }
318  }
319  }
320  catch(const exception& error) {}
321 
322  } while (next_permutation(data.begin(), __end1, __end, hit_type::compare));
323 
324  ymin -= STANDARD_DEVIATIONS * STANDARD_DEVIATIONS;
325  }
326 
327  } else {
328 
329  const int number_of_outliers = NDF - 1;
330 
331  try {
332 
333  (*this)(data.begin(), __end);
334 
335  V.set(*this, data.begin(), __end, gridAngle_deg, sigma_ns);
336  Y.set(*this, data.begin(), __end);
337 
338  V.invert();
339 
340  for (int n = 0; n <= number_of_outliers; ++n) {
341 
342  double ymax = 0.0;
343  int k = -1;
344 
345  for (size_t i = 0; i != Y.size(); ++i) {
346 
347  double y = getChi2(Y, V, i);
348 
349  if (y > ymax) {
350  ymax = y;
351  k = i;
352  }
353  }
354 
355  if (ymax < STANDARD_DEVIATIONS * STANDARD_DEVIATIONS) {
356  break;
357  }
358 
359  try {
360 
361  V.update(k, +HIT_OFF);
362 
363  this->update(data.begin(), __end, V);
364 
365  Y.set(*this, data.begin(), __end);
366 
367  tz = *this;
368  NDF -= 1;
369  N -= getCount(data[k]);
370  }
371  catch(const exception& error) {
372 
373  V.update(k, -HIT_OFF);
374 
375  static_cast<JLine1Z&>(*this) = tz;
376 
377  Y.set(*this, data.begin(), __end);
378 
379  break;
380  }
381  }
382 
383  chi2 = getChi2(Y, V);
384  tz = *this;
385  }
386  catch(const exception& error) {}
387  }
388 
389  if (chi2 != numeric_limits<double>::max()) {
390 
391  tz.rotate_back(R);
392 
393  out.push_back(getFit(JHistory(JMUONPREFIT), tz, *dir, getQuality(chi2, N, NDF), NDF));
394  }
395  }
396 
397 
398  if (numberOfPrefits > 0) {
399 
400  JEvt::iterator __end = out.end();
401 
402  if (Qmin > 0) { // sort distinct maxima
403 
404  __end = out.begin(); // output
405 
406  for (JEvt::iterator p1 = out.begin(); p1 != out.end() && (size_t) distance(out.begin(), __end) < numberOfPrefits; ) {
407 
408  JEvt::iterator p2 = p1;
409 
410  for (JEvt::iterator i = p1; i != out.end(); ++i) {
411  if (i->getQ() > p2->getQ()) {
412  p2 = i;
413  }
414  }
415 
416  swap(*p2, *p1);
417 
418  p2 = p1++;
419 
420  sort(p1, out.end(), JPointing(*p2));
421 
422  for (double Q = p2->getQ(); p1 != out.end() && (p1->getQ() >= p2->getQ() - Qmin || p1->getQ() <= Q); Q = (p1++)->getQ()) {}
423 
424  swap(*(__end++), *p2);
425  }
426 
427  } else if (numberOfPrefits < out.size()) { // sort subset
428 
429  advance(__end = out.begin(), numberOfPrefits);
430 
431  partial_sort(out.begin(), __end, out.end(), qualitySorter);
432 
433  } else { // sort all
434 
435  sort(out.begin(), __end, qualitySorter);
436  }
437 
438  // add downward pointing solutions if available but not yet sufficient
439 
440  int nz = numberOfDZMax - count_if(out.begin(), __end, make_predicate(&JFit::getDZ, DZMax, JComparison::le()));
441 
442  if (nz > 0) {
443 
444  JEvt::iterator __p = __end;
445  JEvt::iterator __q = __end = partition(__p, out.end(), make_predicate(&JFit::getDZ, DZMax, JComparison::le()));
446 
447  if (nz < distance(__p, __q)) {
448 
449  advance(__end = __p, nz);
450 
451  partial_sort(__p, __end, __q, qualitySorter);
452 
453  } else {
454 
455  sort(__p, __end, qualitySorter);
456  }
457  }
458 
459  out.erase(__end, out.end());
460 
461  } else {
462 
463  sort(out.begin(), out.end(), qualitySorter);
464  }
465 
466  return out;
467  }
468 
469 
470  /**
471  * Auxiliary data structure for sorting of hits.
472  */
473  static const struct cmz {
474  /**
475  * Sort hits according times corrected for position along z-axis.
476  *
477  * \param first first hit
478  * \param second second hit
479  * \return true if first hit earlier than second hit; else false
480  */
481  template<class T>
482  inline bool operator()(const T& first, const T& second) const
483  {
484  using namespace JPP;
485 
486  return (first .getT() * getSpeedOfLight() - first .getZ() <
487  second.getT() * getSpeedOfLight() - second.getZ());
488  }
489  } cmz;
490 
491 
494  int debug;
495 
496  private:
497  /**
498  * Configure internal buffer(s).
499  */
500  void configure()
501  {
502  using namespace JPP;
503 
506  }
507 
511  };
512 }
513 
514 #endif
Algorithms for hit clustering and sorting.
TPaveText * p1
Auxiliary methods to evaluate Poisson probabilities and chi2.
Reduced data structure for L1 hit.
Linear fit of JFIT::JLine1Z.
Match operator for Cherenkov light from muon with given direction.
Match operator for Cherenkov light from muon in any direction.
Mathematical constants.
General purpose messaging.
#define WARNING(A)
Definition: JMessage.hh:65
Direct access to module in detector data structure.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Router for direct addressing of module data in detector data structure.
Linear fit of straight line parallel to z-axis to set of hits (objects with position and time).
Template definition of linear fit.
Definition: JEstimator.hh:25
Data structure for set of track fit results.
double getDZ() const
Get Z-slope.
Data structure for fit of straight line paralel to z-axis.
Definition: JLine1Z.hh:29
Determination of the co-variance matrix of hits for a track along z-axis (JFIT::JLine1Z).
Definition: JMatrixNZ.hh:30
void set(const JVector3D &pos, T __begin, T __end, const double alpha, const double sigma)
Set co-variance matrix.
Definition: JMatrixNZ.hh:85
Determination of the time residual vector of hits for a track along z-axis (JFIT::JLine1Z).
Definition: JVectorNZ.hh:23
void set(const JLine1Z &track, T __begin, T __end)
Set time residual vector.
Definition: JVectorNZ.hh:68
Direction set covering (part of) solid angle.
Definition: JOmega3D.hh:68
JPosition3D & rotate_back(const JRotation3D &R)
Rotate back.
Definition: JPosition3D.hh:200
Rotation matrix.
Definition: JRotation3D.hh:114
const JClass_t & getReference() const
Get reference to object.
Definition: JReference.hh:38
Auxiliary class to compare fit results with respect to a reference direction (e.g....
Template L0 hit builder.
Definition: JBuildL0.hh:38
Template L2 builder.
Definition: JBuildL2.hh:49
Reduced data structure for L1 hit.
Definition: JHitR1.hh:35
1D match criterion.
Definition: JMatch1D.hh:33
3D match criterion with road width.
Definition: JMatch3B.hh:36
int getModuleID() const
Get module identifier.
static const int JMUONPREFIT
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition: JHead.cc:162
double getQuality(const double chi2, const int NDF)
Get quality of fit.
JFit getFit(const int id, const JMODEL::JString &string)
Get fit parameters of string.
int getNumberOfPMTs(const JModule &module)
Get number of PMTs.
double getChi2(const double P)
Get chi2 corresponding to given probability.
Definition: JFitToolkit.hh:56
size_t getCount(const array_type< T > &buffer, const JCompare_t &compare)
Count number of unique values.
Definition: JVectorize.hh:261
JPredicate< JResult_t T::*, JComparison::eq > make_predicate(JResult_t T::*member, const JResult_t value)
Helper method to create predicate for data member.
Definition: JPredicate.hh:128
Auxiliary classes and methods for mathematical operations.
Definition: JEigen3D.hh:88
static const double PI
Mathematical constants.
const double getSpeedOfLight()
Get speed of light.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Model fits to data.
JFIT::JHistory JHistory
Definition: JHistory.hh:354
bool qualitySorter(const JFit &first, const JFit &second)
Comparison of fit results.
counter_type advance(counter_type &counter, const counter_type value, const counter_type limit=std::numeric_limits< counter_type >::max())
Advance counter.
const int n
Definition: JPolint.hh:786
bool next_permutation(T __begin, T __last, T __end, JComparator_t compare, std::bidirectional_iterator_tag)
Implementation of method next_permutation for bidirectional iterators.
Definition: JPermutation.hh:20
static const struct JTRIGGER::clusterize clusterize
static const struct JTRIGGER::clusterizeWeight clusterizeWeight
static const JModuleCounter getNumberOfModules
Function object to count unique modules.
Definition: JSTDTypes.hh:14
void update(const size_t k, const double value)
Update inverted matrix at given diagonal element.
Definition: JMatrixNS.hh:446
void invert()
Invert matrix according LDU decomposition.
Definition: JMatrixNS.hh:75
Data structure for fit parameters.
int numberOfOutliers
maximum number of outliers
double ctMin
minimal cosine space angle between PMT axes
double DZMax
maximal slope for downward pointing solutions
int factoryLimit
factory limit for combinatorics
double TMaxLocal_ns
time window for local coincidences [ns]
double gridAngle_deg
grid angle for directions [deg]
size_t numberOfDZMax
additional number of downward pointing solutions
Auxiliary data structure for sorting of hits.
Definition: JMuonPrefit.hh:473
bool operator()(const T &first, const T &second) const
Sort hits according times corrected for position along z-axis.
Definition: JMuonPrefit.hh:482
Wrapper class to make pre-fit of muon trajectory.
Definition: JMuonPrefit.hh:78
static const struct JRECONSTRUCTION::JMuonPrefit::cmz cmz
JTRIGGER::JHitR1 hit_type
Definition: JMuonPrefit.hh:80
JEvt operator()(const KM3NETDAQ::JDAQEvent &event)
Fit function.
Definition: JMuonPrefit.hh:131
JEvt operator()(const buffer_type &dataL0, const buffer_type &dataL1)
Fit function.
Definition: JMuonPrefit.hh:187
JMuonPrefit(const JMuonPrefitParameters_t &parameters, const JModuleRouter &router, const JOmega3D &omega, const int debug=0)
Constructor.
Definition: JMuonPrefit.hh:112
JMuonPrefit(const JMuonPrefitParameters_t &parameters, const JModuleRouter &router, const int debug=0)
Constructor.
Definition: JMuonPrefit.hh:92
const JModuleRouter & router
Definition: JMuonPrefit.hh:492
std::vector< hit_type > buffer_type
Definition: JMuonPrefit.hh:81
void configure()
Configure internal buffer(s).
Definition: JMuonPrefit.hh:500
JEstimator< JLine1Z > JEstimator_t
Definition: JMuonPrefit.hh:79
Auxiliary data structure for sorting of hits.
Definition: JHitR1.hh:203
Data structure for L2 parameters.