Commit f33754b5 authored by David Blyth's avatar David Blyth

saving state

parent 99fe1f87
......@@ -174,27 +174,16 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter
};
while (available_obs.size() >= n_validate) {
auto track = new genfit::Track;
// build genfit::Track object and process with Kalman
// build track seed
std::vector<uint64_t> track_obs;
std::sort(available_obs.begin(), available_obs.end(), seed_sort_fn);
auto measurement = measurements[available_obs.back()];
track->insertMeasurement(measurement->clone());
track_obs.push_back(available_obs.back());
available_obs.pop_back();
auto measurement = measurements[track_obs[0]];
TVector3 ip_pos(0, 0, 0);
auto meas_coords = measurement->getRawHitCoords();
std::cout << meas_coords[0] << "\t" << meas_coords[1] << "\t" << meas_coords[2] << std::endl;
TVector3 meas_pos(0.01 * meas_coords[0], 0.01 * meas_coords[1], 0.01 * meas_coords[2]);
track->setStateSeed(ip_pos, meas_pos);
track->addTrackRep(new genfit::RKTrackRep(11));
track->addTrackRep(new genfit::RKTrackRep(-11));
track->addTrackRep(new genfit::RKTrackRep(13));
track->addTrackRep(new genfit::RKTrackRep(-13));
for (auto rep : track->getTrackReps()) fitter->processTrackPartially(track, rep, 0, 0);
track_obs.push_back(available_obs.back());
available_obs.pop_back();
// extend Track
auto kernel_sort_fn = [&](const uint64_t a_id, const uint64_t b_id) {
auto a_coords = measurements[a_id]->getRawHitCoords();
TVector3 a_pos(a_coords[0], a_coords[1], a_coords[2]);
......@@ -202,30 +191,29 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter
TVector3 b_pos(b_coords[0], b_coords[1], b_coords[2]);
return meas_pos.Dot(a_pos) / a_pos.Mag() > meas_pos.Dot(b_pos) / b_pos.Mag();
};
// auto kernel_sort_fn = [&](const uint64_t a_id, const uint64_t b_id) {
// double meas_param = seed_sort_param[track_obs.back()];
// double a_diff = seed_sort_param[a_id] - meas_param;
// double b_diff = seed_sort_param[b_id] - meas_param;
// if (a_diff * b_diff < 0) return a_diff < b_diff;
// return a_diff > b_diff;
//};
// auto kernel_sort_fn = [&](const uint64_t a_id, const uint64_t b_id) {
// double a_len = 0;
// double b_len = 0;
// for (auto rep : all_reps) {
// auto state = track->getPoint(-1)->getFitterInfo(rep)->getFittedState();
// a_len += rep->extrapolateToPlane(state, measurements[a_id]->constructPlane(state));
// b_len += rep->extrapolateToPlane(state, measurements[b_id]->constructPlane(state));
// }
// if (a_len * b_len < 0) return a_len < b_len;
// return a_len > b_len;
//};
std::sort(available_obs.begin(), available_obs.end(), kernel_sort_fn);
auto track = new genfit::Track;
track->setStateSeed(ip_pos, meas_pos);
track->addTrackRep(new genfit::RKTrackRep(11));
track->addTrackRep(new genfit::RKTrackRep(-11));
track->addTrackRep(new genfit::RKTrackRep(13));
track->addTrackRep(new genfit::RKTrackRep(-13));
int n_seed = 4;
for (int i = 1; i < n_seed; i++) track_obs.push_back(available_obs[i - 1]);
for (auto id : track_obs) {
track->insertMeasurement(measurements[id]->clone());
track->getPoint(-1)->setSortingParameter(seed_sort_param[id]);
}
available_obs.erase(available_obs.begin(), available_obs.begin() + n_seed - 1);
track->sort();
for (auto rep : track->getTrackReps()) fitter->processTrackPartially(track, rep, 0, -1);
for (auto iter = available_obs.begin(); iter != available_obs.end();) {
auto measurement = measurements[*iter];
auto coords = measurement->getRawHitCoords();
std::cout << coords[0] << "\t" << coords[1] << "\t" << coords[2] << std::endl;
track->insertMeasurement(measurements[*iter]->clone(), -1);
track->insertMeasurement(measurements[*iter]->clone());
auto track_point = track->getPoint(-1);
track_point->setSortingParameter(seed_sort_param[*iter]);
......@@ -234,8 +222,7 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter
std::vector<genfit::AbsTrackRep *> bad_reps;
if (track->getNumPoints() >= n_validate) {
for (auto rep : all_reps) {
track->sort();
fitter->processTrackPartially(track, rep, 0, -1);
fitter->processTrackPartially(track, rep, -2, -1);
auto fi = static_cast<genfit::KalmanFitterInfo *>(track_point->getKalmanFitterInfo(rep));
if (!fi) {
bad_reps.push_back(rep);
......@@ -244,7 +231,7 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter
auto update = fi->getUpdate(1);
auto chi2_inc = update->getChiSquareIncrement() / update->getNdf();
std::cout << "chi2_inc: " << chi2_inc << std::endl;
if (chi2_inc > 100) {
if (chi2_inc > 10) {
bad_reps.push_back(rep);
continue;
}
......@@ -255,11 +242,10 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter
//}
}
if (bad_reps.size() == all_reps.size()) {
// track->deletePoint(-1);
// iter++;
// std::cout << "skipping hit" << std::endl;
// continue;
break;
track->deletePoint(-1);
iter++;
std::cout << "skipping hit" << std::endl;
continue;
}
}
......@@ -274,9 +260,9 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter
std::cout << "track_obs.size(): " << track_obs.size() << std::endl;
while (track_obs.size() >= n_validate) {
auto rep = track->getCardinalRep();
// fitter->processTrackWithRep(track, rep);
fitter->processTrackWithRep(track, rep);
if (true /*track->getFitStatus(rep)->isFitConverged()*/) {
if (track->getFitStatus(rep)->isFitConverged()) {
auto eic_track = new eic::Track;
try {
std::cout << "successful fit" << std::endl;
......@@ -325,25 +311,25 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter
fitter->processTrackPartially(track, rep, 1, -1);
double max_chi2_inc = 0;
int max_chi2_inc_i = 1;
for (int i = 1; i < track_obs.size(); i++) {
int max_chi2_inc_i = 0;
for (int i = 0; i < track_obs.size(); i++) {
auto fi =
static_cast<genfit::KalmanFitterInfo *>(track->getPoint(i)->getKalmanFitterInfo(rep));
if (!fi) {
max_chi2_inc_i = i;
break;
}
// auto update = fi->getUpdate(1);
// if (!update) {
// max_chi2_inc_i = i;
// break;
//}
// auto chi2_inc = update->getChiSquareIncrement();
auto chi2_inc = fi->getSmoothedChi2();
if (chi2_inc > max_chi2_inc) {
max_chi2_inc = chi2_inc;
auto update = fi->getUpdate(1);
if (!update) {
max_chi2_inc_i = i;
break;
}
auto chi2_inc = update->getChiSquareIncrement();
// auto chi2_inc = fi->getSmoothedChi2();
// if (chi2_inc > max_chi2_inc) {
// max_chi2_inc = chi2_inc;
// max_chi2_inc_i = i;
//}
}
if (max_chi2_inc_i == 0) {
......@@ -378,12 +364,12 @@ void printUsage() {
}
int main(int argc, char **argv) {
// class NullBuffer : public std::streambuf {
// public:
// int overflow(int c) { return c; }
//};
// NullBuffer nullBuffer;
// genfit::errorOut.rdbuf(&nullBuffer);
class NullBuffer : public std::streambuf {
public:
int overflow(int c) { return c; }
};
NullBuffer nullBuffer;
genfit::errorOut.rdbuf(&nullBuffer);
genfit::MaterialEffects::getInstance()->setNoEffects();
genfit::MaterialEffects::getInstance()->setEnergyLossBrems(false);
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment