Commit a034170d authored by David Blyth's avatar David Blyth

Saving state

parent c31c2cc8
......@@ -156,14 +156,16 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter
if (measurement) {
available_obs.push_back(id);
measurements[id] = measurement;
double t_sum = 0;
double weight_sum = 0;
for (auto pos : obs->pos()) {
double weight = pos.weightmod() + 1;
t_sum += pos.mean().t() * weight;
weight_sum += weight;
}
seed_sort_param[id] = t_sum / weight_sum;
// double t_sum = 0;
// double weight_sum = 0;
// for (auto pos : obs->pos()) {
// double weight = pos.weightmod() + 1;
// t_sum += pos.mean().t() * weight;
// weight_sum += weight;
//}
// seed_sort_param[id] = t_sum / weight_sum;
seed_sort_param[id] =
static_cast<eic::SimHit *>(event->GetEntry(obs->source(0)))->globalprepos().t();
// seed_sort_param[id] = measurement->getRawHitCoords().Norm2Sqr();
}
}
......@@ -182,7 +184,6 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter
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]);
auto kernel_sort_fn = [&](const uint64_t a_id, const uint64_t b_id) {
auto a_coords = measurements[a_id]->getRawHitCoords();
......@@ -196,13 +197,13 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter
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;
int n_seed = 5;
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());
meas_coords = measurements[id]->getRawHitCoords();
std::cout << meas_coords[0] << "\t" << meas_coords[1] << "\t" << meas_coords[2] << std::endl;
track->getPoint(-1)->setSortingParameter(seed_sort_param[id]);
}
available_obs.erase(available_obs.begin(), available_obs.begin() + n_seed - 1);
......@@ -210,48 +211,86 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter
for (auto rep : track->getTrackReps()) fitter->processTrackPartially(track, rep, 0, -1);
for (auto iter = available_obs.begin(); iter != available_obs.end();) {
auto tmp_track = new genfit::Track;
tmp_track->setStateSeed(ip_pos, meas_pos);
tmp_track->addTrackRep(new genfit::RKTrackRep(11));
tmp_track->addTrackRep(new genfit::RKTrackRep(13));
for (auto id : track_obs) {
tmp_track->insertMeasurement(measurements[id]->clone());
tmp_track->getPoint(-1)->setSortingParameter(seed_sort_param[id]);
}
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());
auto track_point = track->getPoint(-1);
tmp_track->insertMeasurement(measurements[*iter]->clone());
auto track_point = tmp_track->getPoint(-1);
track_point->setSortingParameter(seed_sort_param[*iter]);
tmp_track->sort();
// process and check for goodness
std::vector<genfit::AbsTrackRep *> all_reps(track->getTrackReps());
std::vector<genfit::AbsTrackRep *> all_reps(tmp_track->getTrackReps());
std::vector<genfit::AbsTrackRep *> bad_reps;
if (track->getNumPoints() >= n_validate) {
if (tmp_track->getNumPoints() >= n_validate) {
for (auto rep : all_reps) {
fitter->processTrackPartially(track, rep, -2, -1);
auto fi = static_cast<genfit::KalmanFitterInfo *>(track_point->getKalmanFitterInfo(rep));
if (!fi) {
bad_reps.push_back(rep);
continue;
}
auto update = fi->getUpdate(1);
auto chi2_inc = update->getChiSquareIncrement() / update->getNdf();
std::cout << "chi2_inc: " << chi2_inc << std::endl;
if (chi2_inc > 10) {
bad_reps.push_back(rep);
continue;
fitter->processTrackPartially(tmp_track, rep, 0, -1);
tmp_track->reverseTrack();
fitter->processTrackPartially(tmp_track, rep, 0, -1);
tmp_track->reverseTrack();
fitter->processTrackPartially(tmp_track, rep, 0, -1);
// auto fi = static_cast<genfit::KalmanFitterInfo
// *>(track_point->getKalmanFitterInfo(rep)); if (!fi) {
// bad_reps.push_back(rep);
// continue;
//}
// auto update = fi->getUpdate(1);
// auto chi2_inc = update->getChiSquareIncrement() / update->getNdf();
// std::cout << "chi2_inc for pdg code " << rep->getPDG() << ": " << chi2_inc <<
// std::endl; if (chi2_inc > 5) {
// bad_reps.push_back(rep);
// continue;
//}
// try {
// auto pval = fitter->getPVal(tmp_track, rep, 1);
// std::cout << "pval for pdg code " << rep->getPDG() << ": " << pval << std::endl;
// if (fabs(pval - 0.5) < .49) continue;
//} catch (genfit::Exception e) {
// ;
//}
double chi2 = 0;
double ndf = -1. * rep->getDim();
auto points = tmp_track->getPointsWithMeasurement();
for (auto iter = points.begin() + 0; iter != points.end(); iter++) {
auto fi = (*iter)->getKalmanFitterInfo(rep);
if (!fi) continue;
auto update = fi->getForwardUpdate();
chi2 += update->getChiSquareIncrement();
std::cout << "chi2 inc: " << update->getChiSquareIncrement() << std::endl;
ndf += update->getNdf();
}
// auto len = track->getTrackLen(rep, -2, -1);
std::cout << "chi2: " << chi2 << ", ndf: " << ndf << std::endl;
double pval = std::max(0., ROOT::Math::chisquared_cdf_c(chi2, ndf));
std::cout << "pval for pdg code " << rep->getPDG() << ": " << pval << std::endl;
if (fabs(pval - 0.5) < 0.49) continue;
bad_reps.push_back(rep);
// auto len = tmp_track->getTrackLen(rep, -2, -1);
// if (len < 0) {
// bad_reps.push_back(rep);
// continue;
//}
}
if (bad_reps.size() == all_reps.size()) {
track->deletePoint(-1);
iter++;
std::cout << "skipping hit" << std::endl;
continue;
break;
}
}
// clear bad track reps and pop hit from available list
for (auto rep : bad_reps) track->deleteTrackRep(track->getIdForRep(rep));
for (auto rep : bad_reps) tmp_track->deleteTrackRep(tmp_track->getIdForRep(rep));
track_obs.push_back(*iter);
delete track;
track = tmp_track;
iter = available_obs.erase(iter);
}
track->determineCardinalRep();
......@@ -307,8 +346,6 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter
break;
} else {
std::cout << "failed to converge" << std::endl;
track->deleteFitterInfo(0, -1, rep);
fitter->processTrackPartially(track, rep, 1, -1);
double max_chi2_inc = 0;
int max_chi2_inc_i = 0;
......
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