Jpp  18.0.0-rc.2
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  advance(__end = data.begin(), NMaxHits);
220 
221  partial_sort(data.begin(), __end, data.end(), cmz);
222  }
223 
224 
225  // 1D cluster
226 
227  __end = clusterizeWeight(data.begin(), __end, match1D);
228 
229  if (useL0) {
230 
231  buffer_type::iterator p = __end; // begin L0 data
232  buffer_type::iterator q = copy(dataL0.begin(), dataL0.end(), p); // end L0 data
233 
234  for (buffer_type::iterator i = p; i != q; ++i) {
235 
236  if (find_if(data.begin(), __end, make_predicate(&hit_type::getModuleID, i->getModuleID())) == __end) {
237 
238  i->rotate(R);
239 
240  if (match1D.count(*i, data.begin(), __end) != 0) {
241  *p = *i;
242  ++p;
243  }
244  }
245  }
246 
247  __end = clusterize(__end, p, match1D);
248  }
249 
250 
251  if (distance(data.begin(), __end) <= NUMBER_OF_PARAMETERS) {
252  continue;
253  }
254 
255 
256  // 1D fit
257 
258  JLine1Z tz;
259  double chi2 = numeric_limits<double>::max();
260  int NDF = distance(data.begin(), __end) - NUMBER_OF_PARAMETERS;
261  int N = getCount(data.begin(), __end);
262 
263 
264  if (distance(data.begin(), __end) <= factoryLimit) {
265 
266  int number_of_outliers = numberOfOutliers;
267 
268  if (number_of_outliers > NDF - 1) {
269  number_of_outliers = NDF - 1;
270  }
271 
272  double ymin = numeric_limits<double>::max();
273 
274  buffer_type::iterator __end1 = __end;
275 
276  for (int n = 0; n <= number_of_outliers; ++n, --__end1) {
277 
278  sort(data.begin(), __end, hit_type::compare);
279 
280  do {
281  /*
282  if (getNumberOfStrings(router, data.begin(), __end1) < 2) {
283  continue;
284  }
285  */
286  try {
287 
288  (*this)(data.begin(), __end1);
289 
290  V.set(*this, data.begin(), __end1, gridAngle_deg, sigma_ns);
291  Y.set(*this, data.begin(), __end1);
292 
293  V.invert();
294 
295  double y = getChi2(Y, V);
296 
297  if (y <= -(STANDARD_DEVIATIONS * STANDARD_DEVIATIONS)) {
298 
299  WARNING(endl << "chi2(1) " << y << endl);
300 
301  } else {
302 
303  if (y < 0.0) {
304  y = 0.0;
305  }
306 
307  if (y < ymin) {
308  ymin = y;
309  tz = *this;
310  chi2 = ymin;
311  NDF = distance(data.begin(), __end1) - NUMBER_OF_PARAMETERS;
312  N = getCount(data.begin(), __end1);
313  }
314  }
315  }
316  catch(const JException& error) {}
317 
318  } while (next_permutation(data.begin(), __end1, __end, hit_type::compare));
319 
320  ymin -= STANDARD_DEVIATIONS * STANDARD_DEVIATIONS;
321  }
322 
323  } else {
324 
325  const int number_of_outliers = NDF - 1;
326 
327  try {
328 
329  (*this)(data.begin(), __end);
330 
331  V.set(*this, data.begin(), __end, gridAngle_deg, sigma_ns);
332  Y.set(*this, data.begin(), __end);
333 
334  V.invert();
335 
336  for (int n = 0; n <= number_of_outliers; ++n) {
337 
338  double ymax = 0.0;
339  int k = -1;
340 
341  for (size_t i = 0; i != Y.size(); ++i) {
342 
343  double y = getChi2(Y, V, i);
344 
345  if (y > ymax) {
346  ymax = y;
347  k = i;
348  }
349  }
350 
351  if (ymax < STANDARD_DEVIATIONS * STANDARD_DEVIATIONS) {
352  break;
353  }
354 
355  V.update(k, HIT_OFF);
356 
357  this->update(data.begin(), __end, V);
358 
359  Y.set(*this, data.begin(), __end);
360 
361  NDF -= 1;
362  N -= getCount(data[k]);
363  }
364 
365  chi2 = getChi2(Y, V);
366  tz = *this;
367  }
368  catch(const JException& error) {}
369  }
370 
371  if (chi2 != numeric_limits<double>::max()) {
372 
373  tz.rotate_back(R);
374 
375  out.push_back(getFit(JHistory(JMUONPREFIT), tz, *dir, getQuality(chi2, N, NDF), NDF));
376  }
377  }
378 
379  JMuonPrefitParameters_t parameters = static_cast<const JMuonPrefitParameters_t&>(*this);
380 
381  if (parameters.numberOfPrefits > 0) {
382 
383  // apply default sorter
384 
385  JEvt::iterator __end = out.end();
386 
387  if (parameters.numberOfPrefits < out.size()) {
388 
389  advance(__end = out.begin(), parameters.numberOfPrefits);
390 
391  partial_sort(out.begin(), __end, out.end(), qualitySorter);
392 
393  } else {
394 
395  sort(out.begin(), __end, qualitySorter);
396  }
397 
398  // add downward pointing solutions if available but not yet sufficient
399 
400  int nz = parameters.numberOfDZMax - count_if(out.begin(), __end, make_predicate(&JFit::getDZ, parameters.DZMax, JComparison::le()));
401 
402  if (nz > 0) {
403 
404  JEvt::iterator __p = __end;
405  JEvt::iterator __q = __end = partition(__p, out.end(), make_predicate(&JFit::getDZ, parameters.DZMax, JComparison::le()));
406 
407  if (nz < distance(__p, __q)) {
408 
409  advance(__end = __p, nz);
410 
411  partial_sort(__p, __end, __q, qualitySorter);
412 
413  } else {
414 
415  sort(__p, __end, qualitySorter);
416  }
417  }
418 
419  out.erase(__end, out.end());
420 
421  } else {
422 
423  sort(out.begin(), out.end(), qualitySorter);
424  }
425 
426  return out;
427  }
428 
429 
430  /**
431  * Auxiliary data structure for sorting of hits.
432  */
433  static const struct cmz {
434  /**
435  * Sort hits according times corrected for position along z-axis.
436  *
437  * \param first first hit
438  * \param second second hit
439  * \return true if first hit earlier than second hit; else false
440  */
441  template<class T>
442  inline bool operator()(const T& first, const T& second) const
443  {
444  using namespace JPP;
445 
446  return (first .getT() * getSpeedOfLight() - first .getZ() <
447  second.getT() * getSpeedOfLight() - second.getZ());
448  }
449  } cmz;
450 
451 
454  int debug;
455 
456  private:
457  /**
458  * Configure internal buffer(s).
459  */
460  void configure()
461  {
462  using namespace JPP;
463 
466  }
467 
471  };
472 }
473 
474 #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:23
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:442
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:433
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.
Template definition of linear fit.
Definition: JEstimator.hh:25
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:457
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:697
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:80
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:36
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:452
Reduced data structure for L1 hit.
Definition: JHitR1.hh:31
void configure()
Configure internal buffer(s).
Definition: JMuonPrefit.hh:460
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
then echo WARNING
Definition: JTuneHV.sh:91
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