Fit function.
220 {
223
224 const double STANDARD_DEVIATIONS = 3.0;
226
228
231
232 data.reserve(dataL0.size() +
233 dataL1.size());
234
236
238
239 for (JOmega3D_t::const_iterator dir =
omega.begin(); dir !=
omega.end(); ++dir) {
240
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
252
254
255 buffer_type::iterator __p =
data.begin();
256
258
259 partial_sort(
data.begin(), __p, __end,
cmz);
260
261 __end = __p;
262 }
263
264
265
266
268
270
271 buffer_type::iterator p = __end;
272 buffer_type::iterator q =
copy(dataL0.begin(), dataL0.end(), p);
273
274 for (buffer_type::iterator i = p; i != q; ++i) {
275
277
278 i->rotate(R);
279
280 if (match1D.count(*i,
data.begin(), __end) != 0) {
281 *p = *i;
282 ++p;
283 }
284 }
285 }
286
288 }
289
290
291 if (
distance(
data.begin(), __end) <= NUMBER_OF_PARAMETERS) {
292 continue;
293 }
294
295
296
297
299 double chi2 = numeric_limits<double>::max();
300 int NDF =
distance(
data.begin(), __end) - NUMBER_OF_PARAMETERS;
302
303
305
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
319
320 do {
321
322
323
324
325
326 try {
327
328 (*this)(
data.begin(), __end1);
329
332
334
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) {
345 }
346
347 if (y < ymin) {
349 tz = *this;
350 chi2 = ymin;
351 NDF =
distance(
data.begin(), __end1) - NUMBER_OF_PARAMETERS;
353 }
354 }
355 }
356 catch(const exception& error) {}
357
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
373
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
384
385 if (y > ymax) {
387 k = i;
388 }
389 }
390
391 if (ymax < STANDARD_DEVIATIONS * STANDARD_DEVIATIONS) {
392 break;
393 }
394
395 try {
396
398
399 this->update(
data.begin(), __end,
V);
400
402
403 tz = *this;
404 NDF -= 1;
406 }
407 catch(const exception& error) {
408
410
411 static_cast<JLine1Z&
>(*this) = tz;
412
414
415 break;
416 }
417 }
418
420 tz = *this;
421 }
422 catch(const exception& error) {}
423 }
424
425 if (chi2 != numeric_limits<double>::max()) {
426
428
430
431
432
435 }
436 }
437
438
440
441 JEvt::iterator __end = out.end();
442
444
445 __end = out.begin();
446
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
458
460
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
469
471
473
474 } else {
475
477 }
478
479
480
482
483 if (nz > 0) {
484
485 JEvt::iterator __p = __end;
487
489
491
493
494 } else {
495
497 }
498 }
499
500 out.erase(__end, out.end());
501
502 } else {
503
505 }
506
507 return out;
508 }
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
double getDZ() const
Get Z-slope.
Data structure for fit of straight line paralel to z-axis.
void set(const JVector3D &pos, T __begin, T __end, const double alpha, const double sigma)
Set co-variance matrix.
void set(const JLine1Z &track, T __begin, T __end)
Set time residual vector.
JPosition3D & rotate_back(const JRotation3D &R)
Rotate back.
Auxiliary class to compare fit results with respect to a reference direction (e.g....
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.
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.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
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.
JHitIterator_t clusterize(JHitIterator_t __begin, JHitIterator_t __end, const JMatch_t &match, const int Nmin=1)
Partition data according given binary match operator.
void update(const size_t k, const double value)
Update inverted matrix at given diagonal element.
void invert()
Invert matrix according LDU decomposition.
size_t numberOfPrefits
number of prefits
int numberOfOutliers
maximum number of outliers
double DZMax
maximal slope for downward pointing solutions
double Qmin
minimal quality step
int factoryLimit
factory limit for combinatorics
size_t numberOfDZMax
additional number of downward pointing solutions
int NMaxHits
maximal number of hits
double sigma_ns
time resolution [ns]
static const struct JRECONSTRUCTION::JMuonPrefit::cmz cmz