Jpp test-rotations-new
the software that should make you happy
Loading...
Searching...
No Matches
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"
17#include "JTrigger/JMatch1D.hh"
18#include "JTrigger/JMatch3B.hh"
19
20#include "JFit/JLine1Z.hh"
22#include "JFit/JMatrixNZ.hh"
23#include "JFit/JVectorNZ.hh"
24#include "JFit/JFitToolkit.hh"
25
29
30#include "JMath/JConstants.hh"
32
33#include "JLang/JPredicate.hh"
34
37
39
42
43#include "Jeep/JMessage.hh"
44
45
46/**
47 * \author mdejong, gmaggi, azegarelli
48 */
49
50namespace JRECONSTRUCTION {}
51namespace JPP { using namespace JRECONSTRUCTION; }
52
53namespace JRECONSTRUCTION {
54
57 using JFIT::JLine1Z;
58 using JFIT::JEstimator;
59 using JFIT::JMatrixNZ;
60 using JFIT::JVectorNZ;
65
66
67 /**
68 * Wrapper class to make pre-fit of muon trajectory.
69 *
70 * The JMuonPrefit fit is used to generate start values for subsequent fits (usually JMuonSimplex and JMuonGandalf).\n
71 * 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
72 * The directions are spaced by the parameters JMuonPrefitParameters_t::gridAngle_deg.\n
73 * This angle corresponds to the envisaged angular accuracy of the result.\n
74 * The probability that one of the results is less than this angle away from the correct value,
75 * multiple start values should be considered (JMuonPrefitParameters_t::numberOfPrefits).\n
76 * Note that the CPU time scales with the inverse of the square of this angle.\n
77 * The chi-squared is based on the time residuals.\n
78 */
79 struct JMuonPrefit :
81 public JEstimator<JLine1Z>
82 {
86
87 using JEstimator_t::operator();
88
89 /**
90 * Input data type.
91 */
92 struct input_type :
93 public JDAQEventHeader
94 {
95 /**
96 * Default constructor.
97 */
99 {}
100
101
102 /**
103 * Constructor.
104 *
105 * \param header header
106 * \param coverage coverage
107 */
109 const coverage_type& coverage) :
110 JDAQEventHeader(header),
112 {}
113
117 };
118
119
120 /**
121 * Constructor
122 *
123 * \param parameters parameters
124 * \param debug debug
125 */
127 const int debug = 0) :
128 JMuonPrefitParameters_t(parameters),
129 omega (parameters.gridAngle_deg * JMATH::PI/180.0),
130 debug (debug)
131 {}
132
133
134 /**
135 * Constructor
136 *
137 * \param parameters parameters
138 * \param omega directions
139 * \param debug debug
140 */
142 const JOmega3D& omega,
143 const int debug = 0) :
144 JMuonPrefitParameters_t(parameters),
145 omega (omega),
146 debug (debug)
147 {}
148
149
150 /**
151 * Get input data.
152 *
153 * \param router module router
154 * \param event event
155 * \param coverage coverage
156 * \return input data
157 */
159 const JDAQEvent& event,
160 const coverage_type& coverage) const
161
162 {
163 using namespace std;
164 using namespace JTRIGGER;
165 using namespace KM3NETDAQ;
166
167 const JBuildL0<hit_type> buildL0;
169
170 input_type input(event.getDAQEventHeader(), coverage);
171
172 buffer_type& dataL0 = input.dataL0;
173 buffer_type& dataL1 = input.dataL1;
174
175 buildL2(event, router, !useL0, back_inserter(dataL1));
176
177 // 3D cluster of unique optical modules
178
180
181 sort(dataL1.begin(), dataL1.end(), hit_type::compare);
182
183 buffer_type::iterator __end = dataL1.end();
184
185 __end = unique(dataL1.begin(), __end, equal_to<JDAQModuleIdentifier>());
186
187 __end = clusterizeWeight(dataL1.begin(), __end, match3B);
188
189 dataL1.erase(__end, dataL1.end());
190
191
192 if (useL0) {
193
194 buildL0(event, router, true, back_inserter(dataL0));
195
196 __end = dataL0.end();
197
198 for (buffer_type::iterator i = dataL0.begin(); i != __end; ) {
199
200 if (match3B.count(*i, dataL1.begin(), dataL1.end()) != 0)
201 ++i;
202 else
203 swap(*i, *--__end);
204 }
205
206 dataL0.erase(__end, dataL0.end());
207 }
208
209 return input;
210 }
211
212
213 /**
214 * Fit function.
215 *
216 * \param input input data
217 * \return fit results
218 */
220 {
221 using namespace std;
222 using namespace JPP;
223
224 const double STANDARD_DEVIATIONS = 3.0; // [unit]
225 const double HIT_OFF = 1.0e3 * sigma_ns * sigma_ns; // [ns^2]
226
228
229 const buffer_type& dataL0 = input.dataL0;
230 const buffer_type& dataL1 = input.dataL1;
231
232 data.reserve(dataL0.size() +
233 dataL1.size());
234
235 JEvent event(JMUONPREFIT);
236
237 JEvt out;
238
239 for (JOmega3D_t::const_iterator dir = omega.begin(); dir != omega.end(); ++dir) {
240
241 const JRotation3D R(*dir);
242
243
244 buffer_type::iterator __end = copy(dataL1.begin(), dataL1.end(), data.begin());
245
246 for (buffer_type::iterator i = data.begin(); i != __end; ++i) {
247 i->rotate(R);
248 }
249
250
251 // reduce data
252
253 if (distance(data.begin(), __end) > NMaxHits) {
254
255 buffer_type::iterator __p = data.begin();
256
257 advance(__p, NMaxHits);
258
259 partial_sort(data.begin(), __p, __end, cmz);
260
261 __end = __p;
262 }
263
264
265 // 1D cluster
266
267 __end = clusterizeWeight(data.begin(), __end, match1D);
268
269 if (useL0) {
270
271 buffer_type::iterator p = __end; // begin L0 data
272 buffer_type::iterator q = copy(dataL0.begin(), dataL0.end(), p); // end L0 data
273
274 for (buffer_type::iterator i = p; i != q; ++i) {
275
276 if (find_if(data.begin(), __end, make_predicate(&hit_type::getModuleID, i->getModuleID())) == __end) {
277
278 i->rotate(R);
279
280 if (match1D.count(*i, data.begin(), __end) != 0) {
281 *p = *i;
282 ++p;
283 }
284 }
285 }
286
287 __end = clusterize(__end, p, match1D);
288 }
289
290
291 if (distance(data.begin(), __end) <= NUMBER_OF_PARAMETERS) {
292 continue;
293 }
294
295
296 // 1D fit
297
298 JLine1Z tz;
299 double chi2 = numeric_limits<double>::max();
300 int NDF = distance(data.begin(), __end) - NUMBER_OF_PARAMETERS;
301 int N = getCount(data.begin(), __end);
302
303
304 if (distance(data.begin(), __end) <= factoryLimit) {
305
306 int number_of_outliers = numberOfOutliers;
307
308 if (number_of_outliers > NDF - 1) {
309 number_of_outliers = NDF - 1;
310 }
311
312 double ymin = numeric_limits<double>::max();
313
314 buffer_type::iterator __end1 = __end;
315
316 for (int n = 0; n <= number_of_outliers; ++n, --__end1) {
317
318 sort(data.begin(), __end, hit_type::compare);
319
320 do {
321 /*
322 if (getNumberOfStrings(router, data.begin(), __end1) < 2) {
323 continue;
324 }
325 */
326 try {
327
328 (*this)(data.begin(), __end1);
329
330 V.set(*this, data.begin(), __end1, gridAngle_deg, sigma_ns);
331 Y.set(*this, data.begin(), __end1);
332
333 V.invert();
334
335 double y = getChi2(Y, V);
336
337 if (y <= -(STANDARD_DEVIATIONS * STANDARD_DEVIATIONS)) {
338
339 WARNING(endl << "chi2(1) " << y << endl);
340
341 } else {
342
343 if (y < 0.0) {
344 y = 0.0;
345 }
346
347 if (y < ymin) {
348 ymin = y;
349 tz = *this;
350 chi2 = ymin;
351 NDF = distance(data.begin(), __end1) - NUMBER_OF_PARAMETERS;
352 N = getCount(data.begin(), __end1);
353 }
354 }
355 }
356 catch(const exception& error) {}
357
358 } while (next_permutation(data.begin(), __end1, __end, hit_type::compare));
359
360 ymin -= STANDARD_DEVIATIONS * STANDARD_DEVIATIONS;
361 }
362
363 } else {
364
365 const int number_of_outliers = NDF - 1;
366
367 try {
368
369 (*this)(data.begin(), __end);
370
371 V.set(*this, data.begin(), __end, gridAngle_deg, sigma_ns);
372 Y.set(*this, data.begin(), __end);
373
374 V.invert();
375
376 for (int n = 0; n <= number_of_outliers; ++n) {
377
378 double ymax = 0.0;
379 int k = -1;
380
381 for (size_t i = 0; i != Y.size(); ++i) {
382
383 double y = getChi2(Y, V, i);
384
385 if (y > ymax) {
386 ymax = y;
387 k = i;
388 }
389 }
390
391 if (ymax < STANDARD_DEVIATIONS * STANDARD_DEVIATIONS) {
392 break;
393 }
394
395 try {
396
397 V.update(k, +HIT_OFF);
398
399 this->update(data.begin(), __end, V);
400
401 Y.set(*this, data.begin(), __end);
402
403 tz = *this;
404 NDF -= 1;
405 N -= getCount(data[k]);
406 }
407 catch(const exception& error) {
408
409 V.update(k, -HIT_OFF);
410
411 static_cast<JLine1Z&>(*this) = tz;
412
413 Y.set(*this, data.begin(), __end);
414
415 break;
416 }
417 }
418
419 chi2 = getChi2(Y, V);
420 tz = *this;
421 }
422 catch(const exception& error) {}
423 }
424
425 if (chi2 != numeric_limits<double>::max()) {
426
427 tz.rotate_back(R);
428
429 out.push_back(getFit(event(), tz, *dir, getQuality(chi2, N, NDF), NDF));
430
431 // set additional values
432
433 out.rbegin()->setW(JPP_COVERAGE_ORIENTATION, input.coverage.orientation);
434 out.rbegin()->setW(JPP_COVERAGE_POSITION, input.coverage.position);
435 }
436 }
437
438
439 if (numberOfPrefits > 0) {
440
441 JEvt::iterator __end = out.end();
442
443 if (Qmin > 0) { // sort distinct maxima
444
445 __end = out.begin(); // output
446
447 for (JEvt::iterator p1 = out.begin(); p1 != out.end() && (size_t) distance(out.begin(), __end) < numberOfPrefits; ) {
448
449 JEvt::iterator p2 = p1;
450
451 for (JEvt::iterator i = p1; i != out.end(); ++i) {
452 if (i->getQ() > p2->getQ()) {
453 p2 = i;
454 }
455 }
456
457 swap(*p2, *p1);
458
459 p2 = p1++;
460
461 sort(p1, out.end(), JPointing(*p2));
462
463 for (double Q = p2->getQ(); p1 != out.end() && (p1->getQ() >= p2->getQ() - Qmin || p1->getQ() <= Q); Q = (p1++)->getQ()) {}
464
465 swap(*(__end++), *p2);
466 }
467
468 } else if (numberOfPrefits < out.size()) { // sort subset
469
470 advance(__end = out.begin(), numberOfPrefits);
471
472 partial_sort(out.begin(), __end, out.end(), qualitySorter);
473
474 } else { // sort all
475
476 sort(out.begin(), __end, qualitySorter);
477 }
478
479 // add downward pointing solutions if available but not yet sufficient
480
481 int nz = numberOfDZMax - count_if(out.begin(), __end, make_predicate(&JFit::getDZ, DZMax, JComparison::le()));
482
483 if (nz > 0) {
484
485 JEvt::iterator __p = __end;
486 JEvt::iterator __q = __end = partition(__p, out.end(), make_predicate(&JFit::getDZ, DZMax, JComparison::le()));
487
488 if (nz < distance(__p, __q)) {
489
490 advance(__end = __p, nz);
491
492 partial_sort(__p, __end, __q, qualitySorter);
493
494 } else {
495
496 sort(__p, __end, qualitySorter);
497 }
498 }
499
500 out.erase(__end, out.end());
501
502 } else {
503
504 sort(out.begin(), out.end(), qualitySorter);
505 }
506
507 return out;
508 }
509
510
511 /**
512 * Auxiliary data structure for sorting of hits.
513 */
514 static const struct cmz {
515 /**
516 * Sort hits according times corrected for position along z-axis.
517 *
518 * \param first first hit
519 * \param second second hit
520 * \return true if first hit earlier than second hit; else false
521 */
522 template<class T>
523 inline bool operator()(const T& first, const T& second) const
524 {
525 using namespace JPP;
526
527 return (first .getT() * getSpeedOfLight() - first .getZ() <
528 second.getT() * getSpeedOfLight() - second.getZ());
529 }
531
532
534 int debug;
535
536 private:
540 };
541}
542
543#endif
Algorithms for hit clustering and sorting.
Coverage of dynamical detector calibration.
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.
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.
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
const JDAQEventHeader & getDAQEventHeader() const
Get DAQ event header.
int getModuleID() const
Get module identifier.
static const int JMUONPREFIT
static const int JPP_COVERAGE_POSITION
coverage of dynamic position calibration from any Jpp application
static const int JPP_COVERAGE_ORIENTATION
coverage of dynamic orientation calibration from any Jpp application
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition JHead.cc:163
double getChi2(const double P)
Get chi2 corresponding to given probability.
size_t getCount(const array_type< T > &buffer, const JCompare_t &compare)
Count number of unique values.
Auxiliary classes and methods for mathematical operations.
Definition JEigen3D.hh:88
const double getSpeedOfLight()
Get speed of light.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Model fits to data.
double getQuality(const double chi2, const int N, const int NDF)
Get quality of fit.
bool qualitySorter(const JFit &first, const JFit &second)
Comparison of fit results.
JFit getFit(const JHistory &history, const JTrack3D &track, const double Q, const int NDF, const double energy=0.0, const int status=SINGLE_STAGE)
Get fit.
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:791
bool next_permutation(T __begin, T __last, T __end, JComparator_t compare, std::bidirectional_iterator_tag)
Implementation of method next_permutation for bidirectional iterators.
Auxiliary classes and methods for triggering.
JHitIterator_t clusterize(JHitIterator_t __begin, JHitIterator_t __end, const JMatch_t &match, const int Nmin=1)
Partition data according given binary match operator.
Definition JAlgorithm.hh:38
JHitIterator_t clusterizeWeight(JHitIterator_t __begin, JHitIterator_t __end, const JMatch_t &match)
Partition data according given binary match operator.
KM3NeT DAQ data structures and auxiliaries.
Definition DataQueue.cc:39
Data structure for coverage of detector by dynamical calibrations.
Definition JCoverage.hh:19
double position
coverage of detector by available position calibration [0,1]
Definition JCoverage.hh:42
double orientation
coverage of detector by available orientation calibration [0,1]
Definition JCoverage.hh:41
Auxiliary class for historical event.
Definition JHistory.hh:36
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
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.
bool operator()(const T &first, const T &second) const
Sort hits according times corrected for position along z-axis.
input_type(const JDAQEventHeader &header, const coverage_type &coverage)
Constructor.
Wrapper class to make pre-fit of muon trajectory.
static const struct JRECONSTRUCTION::JMuonPrefit::cmz cmz
input_type getInput(const JModuleRouter &router, const JDAQEvent &event, const coverage_type &coverage) const
Get input data.
JMuonPrefit(const JMuonPrefitParameters_t &parameters, const int debug=0)
Constructor.
JEstimator< JLine1Z > JEstimator_t
JEvt operator()(const input_type &input)
Fit function.
std::vector< hit_type > buffer_type
JMuonPrefit(const JMuonPrefitParameters_t &parameters, const JOmega3D &omega, const int debug=0)
Constructor.
Auxiliary data structure for sorting of hits.
Definition JHitR1.hh:203
Data structure for L2 parameters.