Jpp  18.6.0-rc.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 JException& 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  V.update(k, HIT_OFF);
360 
361  this->update(data.begin(), __end, V);
362 
363  Y.set(*this, data.begin(), __end);
364 
365  NDF -= 1;
366  N -= getCount(data[k]);
367  }
368 
369  chi2 = getChi2(Y, V);
370  tz = *this;
371  }
372  catch(const JException& error) {}
373  }
374 
375  if (chi2 != numeric_limits<double>::max()) {
376 
377  tz.rotate_back(R);
378 
379  out.push_back(getFit(JHistory(JMUONPREFIT), tz, *dir, getQuality(chi2, N, NDF), NDF));
380  }
381  }
382 
383  JMuonPrefitParameters_t parameters = static_cast<const JMuonPrefitParameters_t&>(*this);
384 
385  if (parameters.numberOfPrefits > 0) {
386 
387  // apply default sorter
388 
389  JEvt::iterator __end = out.end();
390 
391  if (parameters.numberOfPrefits < out.size()) {
392 
393  advance(__end = out.begin(), parameters.numberOfPrefits);
394 
395  partial_sort(out.begin(), __end, out.end(), qualitySorter);
396 
397  } else {
398 
399  sort(out.begin(), __end, qualitySorter);
400  }
401 
402  // add downward pointing solutions if available but not yet sufficient
403 
404  int nz = parameters.numberOfDZMax - count_if(out.begin(), __end, make_predicate(&JFit::getDZ, parameters.DZMax, JComparison::le()));
405 
406  if (nz > 0) {
407 
408  JEvt::iterator __p = __end;
409  JEvt::iterator __q = __end = partition(__p, out.end(), make_predicate(&JFit::getDZ, parameters.DZMax, JComparison::le()));
410 
411  if (nz < distance(__p, __q)) {
412 
413  advance(__end = __p, nz);
414 
415  partial_sort(__p, __end, __q, qualitySorter);
416 
417  } else {
418 
419  sort(__p, __end, qualitySorter);
420  }
421  }
422 
423  out.erase(__end, out.end());
424 
425  } else {
426 
427  sort(out.begin(), out.end(), qualitySorter);
428  }
429 
430  return out;
431  }
432 
433 
434  /**
435  * Auxiliary data structure for sorting of hits.
436  */
437  static const struct cmz {
438  /**
439  * Sort hits according times corrected for position along z-axis.
440  *
441  * \param first first hit
442  * \param second second hit
443  * \return true if first hit earlier than second hit; else false
444  */
445  template<class T>
446  inline bool operator()(const T& first, const T& second) const
447  {
448  using namespace JPP;
449 
450  return (first .getT() * getSpeedOfLight() - first .getZ() <
451  second.getT() * getSpeedOfLight() - second.getZ());
452  }
453  } cmz;
454 
455 
458  int debug;
459 
460  private:
461  /**
462  * Configure internal buffer(s).
463  */
464  void configure()
465  {
466  using namespace JPP;
467 
470  }
471 
475  };
476 }
477 
478 #endif
Auxiliary methods to evaluate Poisson probabilities and chi2.
JMuonPrefit(const JMuonPrefitParameters_t &parameters, const JModuleRouter &router, const int debug=0)
Constructor.
Definition: JMuonPrefit.hh:92
General exception.
Definition: JException.hh:24
#define WARNING(A)
Definition: JMessage.hh:65
1D match criterion.
Definition: JMatch1D.hh:31
then fatal No hydrophone data file $HYDROPHONE_TXT fi sort gr k
Wrapper class to make pre-fit of muon trajectory.
Definition: JMuonPrefit.hh:75
Linear fit of straight line parallel to z-axis to set of hits (objects with position and time)...
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
static struct JTRIGGER::clusterizeWeight clusterizeWeight
int getModuleID() const
Get module identifier.
Match operator for Cherenkov light from muon with given direction.
bool operator()(const T &first, const T &second) const
Sort hits according times corrected for position along z-axis.
Definition: JMuonPrefit.hh:446
Template definition of linear fit.
Definition: JEstimator.hh:25
double getQuality(const double chi2, const int NDF)
Get quality of fit.
void set(const JVector3D &pos, T __begin, T __end, const double alpha, const double sigma)
Set co-variance matrix.
Definition: JMatrixNZ.hh:85
JPosition3D & rotate_back(const JRotation3D &R)
Rotate back.
Definition: JPosition3D.hh:200
JEvt operator()(const KM3NETDAQ::JDAQEvent &event)
Fit function.
Definition: JMuonPrefit.hh:131
Algorithms for hit clustering and sorting.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Auxiliary data structure for sorting of hits.
Definition: JMuonPrefit.hh:437
Router for direct addressing of module data in detector data structure.
Rotation matrix.
Definition: JRotation3D.hh:111
Data structure for fit parameters.
int getNumberOfPMTs(const JModule &module)
Get number of PMTs.
*fatal Wrong number of arguments esac JCookie sh typeset Z DETECTOR typeset Z SOURCE_RUN typeset Z TARGET_RUN set_variable PARAMETERS_FILE $WORKDIR parameters
Definition: diff-Tuna.sh:38
Direction set covering (part of) solid angle.
Definition: JOmega3D.hh:66
Linear fit of JFIT::JLine1Z.
static const JModuleCounter getNumberOfModules
Function object to count unique modules.
void update(const size_t k, const double value)
Update inverted matrix at given diagonal element.
Definition: JMatrixNS.hh:446
static struct JRECONSTRUCTION::JMuonPrefit::cmz cmz
bool qualitySorter(const JFit &first, const JFit &second)
Comparison of fit results.
then echo The file $DIR KM3NeT_00000001_00000000 root already please rename or remove it first
int numberOfOutliers
maximum number of outliers
static struct JTRIGGER::JHitR1::compare compare
const int n
Definition: JPolint.hh:786
static const int JMUONPREFIT
Template L2 builder.
Definition: JBuildL2.hh:45
Determination of the time residual vector of hits for a track along z-axis (JFIT::JLine1Z).
Definition: JVectorNZ.hh:21
Match operator for Cherenkov light from muon in any direction.
Determination of the co-variance matrix of hits for a track along z-axis (JFIT::JLine1Z).
Definition: JMatrixNZ.hh:28
Mathematical constants.
JTRIGGER::JHitR1 hit_type
Definition: JMuonPrefit.hh:80
do set_variable OUTPUT_DIRECTORY $WORKDIR T
double gridAngle_deg
grid angle for directions [deg]
JMuonPrefit(const JMuonPrefitParameters_t &parameters, const JModuleRouter &router, const JOmega3D &omega, const int debug=0)
Constructor.
Definition: JMuonPrefit.hh:112
static const double PI
Mathematical constants.
void set(const JLine1Z &track, T __begin, T __end)
Set time residual vector.
Definition: JVectorNZ.hh:68
General purpose messaging.
counter_type advance(counter_type &counter, const counter_type value, const counter_type limit=std::numeric_limits< counter_type >::max())
Advance counter.
JEstimator< JLine1Z > JEstimator_t
Definition: JMuonPrefit.hh:79
double TMaxLocal_ns
time window for local coincidences [ns]
void invert()
Invert matrix according LDU decomposition.
Definition: JMatrixNS.hh:75
Direct access to module in detector data structure.
Data structure for L2 parameters.
double getDZ() const
Get Z-slope.
int factoryLimit
factory limit for combinatorics
Reduced data structure for L1 hit.
then JCookie sh JDataQuality D $DETECTOR_ID R
Definition: JDataQuality.sh:41
then usage $script< input file >[option[primary[working directory]]] nWhere option can be N
Definition: JMuonPostfit.sh:40
JFIT::JHistory JHistory
Definition: JHistory.hh:354
const JClass_t & getReference() const
Get reference to object.
Definition: JReference.hh:38
const double getSpeedOfLight()
Get speed of light.
Data structure for set of track fit results.
JFit getFit(const int id, const JMODEL::JString &string)
Get fit parameters of string.
then if[[!-f $DETECTOR]] then JDetector sh $DETECTOR fi cat $WORKDIR trigger_parameters txt<< EOFtrigger3DMuon.enabled=1;trigger3DMuon.numberOfHits=5;trigger3DMuon.gridAngle_deg=1;ctMin=0.0;TMaxLocal_ns=15.0;EOF set_variable TRIGGEREFFICIENCY_TRIGGERED_EVENTS_ONLY INPUT_FILES=() for((i=1;$i<=$NUMBER_OF_RUNS;++i));do JSirene.sh $DETECTOR $JPP_DATA/genhen.km3net_wpd_V2_0.evt.gz $WORKDIR/sirene_ ${i}.root JTriggerEfficiency.sh $DETECTOR $DETECTOR $WORKDIR/sirene_ ${i}.root $WORKDIR/trigger_efficiency_ ${i}.root $WORKDIR/trigger_parameters.txt $JPP_DATA/PMT_parameters.txt INPUT_FILES+=($WORKDIR/trigger_efficiency_ ${i}.root) done for ANGLE_DEG in $ANGLES_DEG[*];do set_variable SIGMA_NS 3.0 set_variable OUTLIERS 3 set_variable OUTPUT_FILE $WORKDIR/matrix\[${ANGLE_DEG}\deg\].root $JPP_DIR/examples/JReconstruction-f"$INPUT_FILES[*]"-o $OUTPUT_FILE-S ${SIGMA_NS}-A ${ANGLE_DEG}-O ${OUTLIERS}-d ${DEBUG}--!fiif[[$OPTION=="plot"]];then if((0));then for H1 in h0 h1;do JPlot1D-f"$WORKDIR/matrix["${^ANGLES_DEG}" deg].root:${H1}"-y"1 2e3"-Y-L TR-T""-\^"number of events [a.u.]"-> o chi2
Definition: JMatrixNZ.sh:106
Data structure for fit of straight line paralel to z-axis.
Definition: JLine1Z.hh:27
int getCount(const T &hit)
Get hit count.
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition: JHead.cc:162
const JModuleRouter & router
Definition: JMuonPrefit.hh:456
Reduced data structure for L1 hit.
Definition: JHitR1.hh:31
void configure()
Configure internal buffer(s).
Definition: JMuonPrefit.hh:464
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
JEvt operator()(const buffer_type &dataL0, const buffer_type &dataL1)
Fit function.
Definition: JMuonPrefit.hh:187
Template L0 hit builder.
Definition: JBuildL0.hh:35
double getChi2(const double P)
Get chi2 corresponding to given probability.
Definition: JFitToolkit.hh:56
std::vector< hit_type > buffer_type
Definition: JMuonPrefit.hh:81
double ctMin
minimal cosine space angle between PMT axes
static struct JTRIGGER::clusterize clusterize
3D match criterion with road width.
Definition: JMatch3B.hh:34