Fit function.
189 {
192
193 const double STANDARD_DEVIATIONS = 3.0;
195
197
198 data.reserve(dataL0.size() +
199 dataL1.size());
200
202
203 for (JOmega3D_t::const_iterator dir =
omega.begin(); dir !=
omega.end(); ++dir) {
204
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
216
218
219 buffer_type::iterator __p =
data.begin();
220
222
223 partial_sort(
data.begin(), __p, __end,
cmz);
224
225 __end = __p;
226 }
227
228
229
230
232
234
235 buffer_type::iterator p = __end;
236 buffer_type::iterator q =
copy(dataL0.begin(), dataL0.end(), p);
237
238 for (buffer_type::iterator i = p; i != q; ++i) {
239
241
242 i->rotate(R);
243
244 if (match1D.count(*i,
data.begin(), __end) != 0) {
245 *p = *i;
246 ++p;
247 }
248 }
249 }
250
252 }
253
254
255 if (
distance(
data.begin(), __end) <= NUMBER_OF_PARAMETERS) {
256 continue;
257 }
258
259
260
261
263 double chi2 = numeric_limits<double>::max();
264 int NDF =
distance(
data.begin(), __end) - NUMBER_OF_PARAMETERS;
266
267
269
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
283
284 do {
285
286
287
288
289
290 try {
291
292 (*this)(
data.begin(), __end1);
293
296
298
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) {
309 }
310
311 if (y < ymin) {
313 tz = *this;
314 chi2 = ymin;
315 NDF =
distance(
data.begin(), __end1) - NUMBER_OF_PARAMETERS;
317 }
318 }
319 }
320 catch(const exception& error) {}
321
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
337
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
348
349 if (y > ymax) {
351 k = i;
352 }
353 }
354
355 if (ymax < STANDARD_DEVIATIONS * STANDARD_DEVIATIONS) {
356 break;
357 }
358
359 try {
360
362
363 this->update(
data.begin(), __end,
V);
364
366
367 tz = *this;
368 NDF -= 1;
370 }
371 catch(const exception& error) {
372
374
375 static_cast<JLine1Z&
>(*this) = tz;
376
378
379 break;
380 }
381 }
382
384 tz = *this;
385 }
386 catch(const exception& error) {}
387 }
388
389 if (chi2 != numeric_limits<double>::max()) {
390
392
394 }
395 }
396
397
399
400 JEvt::iterator __end = out.end();
401
403
404 __end = out.begin();
405
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
417
419
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
428
430
432
433 } else {
434
436 }
437
438
439
441
442 if (nz > 0) {
443
444 JEvt::iterator __p = __end;
446
448
450
452
453 } else {
454
456 }
457 }
458
459 out.erase(__end, out.end());
460
461 } else {
462
464 }
465
466 return out;
467 }
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
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.
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.
static const struct JTRIGGER::clusterize clusterize
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