Commit 3d34b079 authored by David Blyth's avatar David Blyth

Pretty functional with provided event kernel and with backup kernel

parent 867572ca
......@@ -170,9 +170,9 @@ double forwardPVal(genfit::Track *track, genfit::AbsTrackRep *rep) {
double ndf = -1. * rep->getDim();
for (auto point : track->getPointsWithMeasurement()) {
auto fi = point->getKalmanFitterInfo(rep);
if (!fi) continue;
if (!fi) return 0;
auto update = fi->getForwardUpdate();
if (!update) continue;
if (!update) return 0;
chi2 += update->getChiSquareIncrement();
ndf += update->getNdf();
}
......@@ -233,7 +233,7 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter
// build track seed
std::vector<uint64_t> track_obs;
std::sort(available_obs.begin(), available_obs.end(), seed_sort_fn);
if (seed_sort_param[available_obs.back()] > 1) break;
// if (seed_sort_param[available_obs.back()] > 1) break;
track_obs.push_back(available_obs.back());
available_obs.pop_back();
auto current_id = track_obs.back();
......@@ -261,10 +261,6 @@ 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();
};
if (kernel_from_event)
std::sort(available_obs.begin(), available_obs.end(), event_kernel_sort_fn);
else
std::sort(available_obs.begin(), available_obs.end(), backup_kernel_sort_fn);
auto track = new genfit::Track;
track->setStateSeed(ip_pos, meas_pos);
......@@ -272,21 +268,29 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter
track->getPoint(0)->setSortingParameter(0.0);
track->addTrackRep(new genfit::RKTrackRep(11));
track->addTrackRep(new genfit::RKTrackRep(13));
int n_seed = 5;
std::map<genfit::TrackPoint *, uint64_t> id_lookup;
int n_seed = 4;
for (int i = 0; i < n_seed - 1; i++) {
track_obs.push_back(available_obs.back());
if (kernel_from_event)
std::sort(available_obs.begin(), available_obs.end(), event_kernel_sort_fn);
else
std::sort(available_obs.begin(), available_obs.end(), backup_kernel_sort_fn);
auto id = available_obs.back();
track_obs.push_back(id);
available_obs.pop_back();
current_id = id;
meas_coords = measurements[id]->getRawHitCoords();
meas_pos = TVector3(0.01 * meas_coords[0], 0.01 * meas_coords[1], 0.01 * meas_coords[2]);
}
std::map<genfit::TrackPoint *, uint64_t> id_lookup;
for (auto id : track_obs) {
track->insertMeasurement(measurements[id]->clone());
meas_coords = measurements[id]->getRawHitCoords();
auto track_point = track->getPoint(-1);
track_point->setSortingParameter(seed_sort_param[id]);
id_lookup[track_point] = id;
}
track->sort();
fitter->setMaxIterations(16);
fitter->setMaxIterations(4);
try {
fitter->processTrack(track);
} catch (genfit::Exception e) {
......@@ -302,13 +306,26 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter
// extend track
while (do_extend && available_obs.size() > 0) {
meas_coords = track->getPoint(-1)->getRawMeasurement()->getRawHitCoords();
meas_pos = TVector3(0.01 * meas_coords[0], 0.01 * meas_coords[1], 0.01 * meas_coords[2]);
current_id = id_lookup[track->getPoint(-1)];
if (kernel_from_event)
std::sort(available_obs.begin(), available_obs.end(), event_kernel_sort_fn);
else
else {
meas_coords = track->getPoint(-1)->getRawMeasurement()->getRawHitCoords();
meas_pos = TVector3(0.01 * meas_coords[0], 0.01 * meas_coords[1], 0.01 * meas_coords[2]);
std::sort(available_obs.begin(), available_obs.end(), backup_kernel_sort_fn);
}
auto id = available_obs.back();
double kernel_val;
if (kernel[current_id].count(id))
kernel_val = kernel[current_id][id];
else
kernel_val = kernel[id][current_id];
// std::cout << kernel_val << std::endl;
if (kernel_from_event && kernel_val < 0.1) {
break;
}
track_obs.push_back(id);
available_obs.pop_back();
auto measurement = measurements[id];
......@@ -317,18 +334,23 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter
auto track_point = track->getPoint(-1);
id_lookup[track_point] = id;
track_point->setSortingParameter(seed_sort_param[id]);
track->sort();
bool order_changed = track->sort();
// if (order_changed) std::cout << "reordered" << std::endl;
std::vector<genfit::AbsTrackRep *> all_reps(track->getTrackReps());
std::vector<genfit::AbsTrackRep *> bad_reps;
for (auto rep : all_reps) {
try {
fitter->processTrackPartially(track, rep, -2, -1);
if (order_changed)
fitter->processTrackWithRep(track, rep);
else
fitter->processTrackPartially(track, rep, -2, -1);
} catch (genfit::Exception e) {
std::cerr << e.what() << std::endl;
}
double pval = forwardPVal(track, rep);
if (fabs(pval - 0.5) < 0.49) continue;
// std::cout << "pval: " << pval << std::endl;
if (pval > 0.001) continue;
bad_reps.push_back(rep);
}
if (bad_reps.size() == all_reps.size()) {
......@@ -344,19 +366,13 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter
}
for (auto rep : bad_reps) track->deleteTrackRep(track->getIdForRep(rep));
// try {
// quickTrackProcess(fitter, track);
//} catch (genfit::Exception e) {
// std::cerr << e.what() << std::endl;
// break;
//}
}
track->determineCardinalRep();
auto rep = track->getCardinalRep();
// fit track
if (track_obs.size() >= n_validate) {
fitter->setMaxIterations(100);
fitter->setMaxIterations(128);
try {
fitter->processTrackWithRep(track, rep);
} catch (genfit::Exception e) {
......@@ -369,6 +385,8 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter
for (unsigned int id = 0; id < track->getNumPoints() - 1; id++) {
try {
auto point = track->getPoint(id);
auto proio_id = id_lookup[point];
if (proio_id > 0) eic_track->add_observation(proio_id);
auto info = static_cast<genfit::KalmanFitterInfo *>(point->getFitterInfo(rep));
auto state = info->getFittedState();
TVector3 pos = rep->getPos(state);
......@@ -399,13 +417,17 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter
std::cerr << e.what() << std::endl;
}
}
auto point = track->getPoint(track->getNumPoints() - 1);
eic_track->add_observation(id_lookup[point]);
auto entryID = event->AddEntry(eic_track, "Tracker");
event->TagEntry(entryID, "Reconstructed");
event->TagEntry(entryID, "Vis");
track_obs.clear();
std::cout << "successful fit" << std::endl;
}
}
} else
std::cout << "fit failure" << std::endl;
// remove largest outlying hit from track
double max_chi2_inc = 0;
......@@ -454,8 +476,10 @@ int main(int argc, char **argv) {
genfit::MaterialEffects::getInstance()->setNoEffects();
genfit::MaterialEffects::getInstance()->setEnergyLossBrems(false);
genfit::KalmanFitter fitter;
fitter.setDeltaPval(0.1);
fitter.setRelChi2Change(1);
int comp = 0;
int n_validate = 8;
int n_validate = 6;
int opt;
while ((opt = getopt(argc, argv, "c:n:h")) != -1) {
......
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