9 #include <phgenfit/Measurement.h>
14 #include <GenFit/AbsFitterInfo.h>
15 #include <GenFit/AbsHMatrix.h>
16 #include <GenFit/AbsMeasurement.h>
17 #include <GenFit/AbsTrackRep.h>
18 #include <GenFit/DetPlane.h>
19 #include <GenFit/Exception.h>
20 #include <GenFit/FitStatus.h>
21 #include <GenFit/KalmanFittedStateOnPlane.h>
22 #include <GenFit/KalmanFitterInfo.h>
23 #include <GenFit/MeasuredStateOnPlane.h>
24 #include <GenFit/MeasurementOnPlane.h>
25 #include <GenFit/SharedPlanePtr.h>
28 #include <GenFit/TrackPoint.h>
30 #include <TMatrixDSymfwd.h>
31 #include <TMatrixDfwd.h>
33 #include <TMatrixTSym.h>
35 #include <TVectorDfwd.h>
45 #define LogDebug(exp) std::cout << "DEBUG: " << __FILE__ << ": " << __LINE__ << ": " << exp << std::endl
46 #define LogError(exp) std::cout << "ERROR: " << __FILE__ << ": " << __LINE__ << ": " << exp << std::endl
47 #define LogWarning(exp) std::cout << "WARNING: " << __FILE__ << ": " << __LINE__ << ": " << exp << std::endl
49 #define WILD_DOUBLE -999999
57 ofstream fout_matrix(
"matrix.txt");
73 TVectorD seedState(6);
74 TMatrixDSym seedCov(6);
166 std::cout <<
"Track has no TrackPoint with fitterInfo! \n";
170 *(static_cast<genfit::KalmanFitterInfo*>(tp->
getFitterInfo(rep))->getBackwardUpdate())));
178 std::cerr <<
"Exception, next track" << std::endl;
179 std::cerr << e.
what();
213 std::cout <<
"Track has no TrackPoint with fitterInfo! \n";
217 *(static_cast<genfit::KalmanFitterInfo*>(tp->
getFitterInfo(rep))->getBackwardUpdate())));
225 std::cerr <<
"Exception, next track" << std::endl;
226 std::cerr << e.
what();
241 double pathlenth = this->
extrapolateToLine(*state, line_point, line_direction, tr_point_id);
254 std::cout << __LINE__ << std::endl;
257 <<
": tr_point_id: " << tr_point_id
258 <<
": direction: " << direction
261 assert(direction == 1 or direction == -1);
274 bool have_tp_with_fit_info =
false;
275 std::unique_ptr<genfit::MeasuredStateOnPlane> kfsop = NULL;
287 if (dynamic_cast<genfit::KalmanFitterInfo*>(tp->
getFitterInfo(rep)))
291 if (static_cast<genfit::KalmanFitterInfo*>(tp->
getFitterInfo(
293 ->getForwardUpdate())
295 have_tp_with_fit_info =
true;
300 ->getForwardUpdate())));
305 if (static_cast<genfit::KalmanFitterInfo*>(tp->
getFitterInfo(
307 ->getBackwardUpdate())
309 have_tp_with_fit_info =
true;
314 ->getBackwardUpdate())));
320 if (!have_tp_with_fit_info)
323 std::cout << __LINE__ <<
": !have_tp_with_fit_info" << std::endl;
329 if (!kfsop)
return pathlenth;
341 std::cerr << e.
what();
355 assert(direction == 1 or direction == -1);
357 double pathlenth = this->
extrapolateToCylinder(*state, radius, line_point, line_direction, tr_point_id, direction);
368 const std::vector<PHGenFit::Measurement*>& measurements,
369 std::map<
double, std::shared_ptr<PHGenFit2::Track> >& incr_chi2s_new_tracks,
370 const int base_tp_idx,
372 const float blowup_factor,
373 const bool use_fitted_state)
const
378 <<
" : base_tp_idx: " << base_tp_idx
379 <<
" : direction: " << direction
380 <<
" : blowup_factor: " << blowup_factor
381 <<
" : use_fitted_state: " << use_fitted_state
385 if (measurements.size() == 0)
return -1;
389 std::shared_ptr<PHGenFit2::Track> new_track = NULL;
403 std::unique_ptr<genfit::MeasuredStateOnPlane> currentState = NULL;
406 if (track->getNumPointsWithMeasurement() > 0)
408 tp_base = track->getPointWithMeasurement(base_tp_idx);
413 std::cout << __LINE__ <<
": "
414 <<
"newFi: " << newFi << std::endl;
420 track->getCovSeed());
439 if (use_fitted_state)
474 if (blowup_factor > 1)
476 currentState->
blowUpCov(blowup_factor,
true, 1e6);
484 <<
": Fitted state not found!"
486 std::cerr << e.
what() << std::endl;
490 track->getCovSeed());
494 std::cout << __LINE__ << std::endl;
506 new_track->addMeasurement(measurement);
509 std::cout << __LINE__ <<
": clusterIDs size: " << new_track->get_cluster_IDs().size() << std::endl;
510 std::cout << __LINE__ <<
": clusterkeyss size: " << new_track->get_cluster_keys().size() << std::endl;
516 std::cout << __LINE__ << std::endl;
523 <<
": track->getPointWithMeasurement(): " << track->getPointWithMeasurement(-1)
533 std::cout << __LINE__ << std::endl;
535 const std::vector<genfit::AbsMeasurement*>& rawMeasurements =
538 plane = rawMeasurements[0]->constructPlane(*currentState);
550 LogWarning(
"Can not extrapolate to measuremnt with cluster_ID and cluster key: ") << measurement->get_cluster_ID()
551 <<
" " << measurement->get_cluster_key()
557 std::cout << __LINE__ << std::endl;
562 std::cout << __LINE__ << std::endl;
564 TVectorD stateVector(state->
getState());
568 std::cout << __LINE__ << std::endl;
576 for (std::vector<genfit::AbsMeasurement*>::const_iterator
it =
577 rawMeasurements.begin();
578 it != rawMeasurements.end(); ++
it)
581 (*it)->constructMeasurementsOnPlane(*state));
587 std::cout << __LINE__ << std::endl;
594 <<
": size of fi's MeasurementsOnPlane: " << measurements_on_plane.size()
597 for (std::vector<genfit::MeasurementOnPlane*>::const_iterator
it =
598 measurements_on_plane.begin();
599 it != measurements_on_plane.end(); ++
it)
604 const TVectorD& measurement(mOnPlane.
getState());
607 const TMatrixDSym& V(mOnPlane.
getCov());
609 TVectorD res(measurement -
H->Hv(stateVector));
612 std::cout << __LINE__ << std::endl;
622 TMatrixDSym covSumInv(
cov);
637 TMatrixD CHt(
H->MHt(
cov));
638 #ifdef _PRINT_MATRIX_
639 std::cout << __LINE__ <<
": V_{k}:" << std::endl;
641 std::cout << __LINE__ <<
": R_{k}^{-1}:" << std::endl;
643 std::cout << __LINE__ <<
": C_{k|k-1}:" << std::endl;
645 std::cout << __LINE__ <<
": C_{k|k-1} H_{k}^{T} :" << std::endl;
647 std::cout << __LINE__ <<
": K_{k} :" << std::endl;
648 TMatrixD Kk(CHt, TMatrixD::kMult, covSumInv);
650 std::cout << __LINE__ <<
": res:" << std::endl;
654 TMatrixD(CHt, TMatrixD::kMult, covSumInv) * res);
658 covSumInv.Similarity(CHt);
662 std::cout << __LINE__ << std::endl;
673 TVectorD resNew(measurement -
H->Hv(stateVector));
676 TMatrixDSym HCHt(
cov);
692 chi2inc += HCHt.Similarity(resNew);
694 ndfInc += measurement.GetNrows();
696 #ifdef _PRINT_MATRIX_
697 std::cout << __LINE__ <<
": V - HCHt:" << std::endl;
699 std::cout << __LINE__ <<
": resNew:" << std::endl;
704 std::cout << __LINE__ <<
": ndfInc: " << ndfInc << std::endl;
705 std::cout << __LINE__ <<
": chi2inc: " << chi2inc << std::endl;
719 incr_chi2s_new_tracks.insert(std::make_pair(chi2inc, new_track));
734 std::cout <<
"Track has no TrackPoint with fitterInfo! \n";
738 *(static_cast<genfit::KalmanFitterInfo*>(tp->
getFitterInfo(rep))->getBackwardUpdate())));
746 std::cerr <<
"Exception, next track" << std::endl;
747 std::cerr << e.
what();
774 double chi2 = std::numeric_limits<double>::quiet_NaN();
788 double ndf = std::numeric_limits<double>::quiet_NaN();
802 double charge = std::numeric_limits<double>::quiet_NaN();
815 if (!tp_base)
return charge;
835 std::cerr <<
"Track::get_charge - Error - obtaining charge failed. Returning NAN as charge." << std::endl;
843 TVector3
mom(0, 0, 0);
854 if (!tp_base)
return mom;
861 if (!kfi)
return mom;