Commit e04cad71 authored by David Blyth's avatar David Blyth

Added more documentation and slightly better parameterization

parent 0afd1faa
...@@ -179,7 +179,7 @@ double forwardPVal(genfit::Track *track, genfit::AbsTrackRep *rep) { ...@@ -179,7 +179,7 @@ double forwardPVal(genfit::Track *track, genfit::AbsTrackRep *rep) {
return std::max(0., ROOT::Math::chisquared_cdf_c(chi2, ndf)); return std::max(0., ROOT::Math::chisquared_cdf_c(chi2, ndf));
} }
void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter *fitter) { void trackEvent(proio::Event *event, const int n_seed, const int n_validate, genfit::KalmanFitter *fitter) {
// build kernel from event // build kernel from event
std::map<uint64_t, std::map<uint64_t, double>> kernel; std::map<uint64_t, std::map<uint64_t, double>> kernel;
bool kernel_from_event = false; bool kernel_from_event = false;
...@@ -234,7 +234,7 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter ...@@ -234,7 +234,7 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter
// build track seed // build track seed
std::vector<uint64_t> track_obs; std::vector<uint64_t> track_obs;
std::sort(available_obs.begin(), available_obs.end(), seed_sort_fn); std::sort(available_obs.begin(), available_obs.end(), seed_sort_fn);
if (seed_sort_param[available_obs.back()] > 1) break; // limit time window for seed hits if (seed_sort_param[available_obs.back()] > 2) break; // enforce time window for seed hits
track_obs.push_back(available_obs.back()); track_obs.push_back(available_obs.back());
available_obs.pop_back(); available_obs.pop_back();
auto current_id = track_obs.back(); auto current_id = track_obs.back();
...@@ -243,6 +243,7 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter ...@@ -243,6 +243,7 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter
auto meas_coords = measurement->getRawHitCoords(); auto meas_coords = measurement->getRawHitCoords();
TVector3 meas_pos(10 * meas_coords[0], 10 * meas_coords[1], 10 * meas_coords[2]); TVector3 meas_pos(10 * meas_coords[0], 10 * meas_coords[1], 10 * meas_coords[2]);
// set up sort function for kernel and backup kernel
auto event_kernel_sort_fn = [&](const uint64_t a_id, const uint64_t b_id) { auto event_kernel_sort_fn = [&](const uint64_t a_id, const uint64_t b_id) {
double kernel_a, kernel_b; double kernel_a, kernel_b;
if (kernel[current_id].count(a_id)) if (kernel[current_id].count(a_id))
...@@ -267,10 +268,12 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter ...@@ -267,10 +268,12 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter
track->setStateSeed(ip_pos, meas_pos); track->setStateSeed(ip_pos, meas_pos);
track->insertMeasurement(makeIPMeasurement()); track->insertMeasurement(makeIPMeasurement());
track->getPoint(0)->setSortingParameter(0.0); track->getPoint(0)->setSortingParameter(0.0);
// add electron hypothesis
track->addTrackRep(new genfit::RKTrackRep(11)); track->addTrackRep(new genfit::RKTrackRep(11));
// add muon hypothesis
track->addTrackRep(new genfit::RKTrackRep(13)); track->addTrackRep(new genfit::RKTrackRep(13));
std::map<genfit::TrackPoint *, uint64_t> id_lookup; std::map<genfit::TrackPoint *, uint64_t> id_lookup;
int n_seed = 4; // add additional seed hits using the kernel to sort by ascending similarity
for (int i = 0; i < n_seed - 1; i++) { for (int i = 0; i < n_seed - 1; i++) {
if (kernel_from_event) if (kernel_from_event)
std::sort(available_obs.begin(), available_obs.end(), event_kernel_sort_fn); std::sort(available_obs.begin(), available_obs.end(), event_kernel_sort_fn);
...@@ -280,6 +283,7 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter ...@@ -280,6 +283,7 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter
track_obs.push_back(id); track_obs.push_back(id);
available_obs.pop_back(); available_obs.pop_back();
// update current hit variables to update the kernel sort functions
current_id = id; current_id = id;
meas_coords = measurements[id]->getRawHitCoords(); meas_coords = measurements[id]->getRawHitCoords();
meas_pos = TVector3(0.01 * meas_coords[0], 0.01 * meas_coords[1], 0.01 * meas_coords[2]); meas_pos = TVector3(0.01 * meas_coords[0], 0.01 * meas_coords[1], 0.01 * meas_coords[2]);
...@@ -307,6 +311,7 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter ...@@ -307,6 +311,7 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter
// extend track // extend track
while (do_extend && available_obs.size() > 0) { while (do_extend && available_obs.size() > 0) {
// sort available hits relative to last hit in track
current_id = id_lookup[track->getPoint(-1)]; current_id = id_lookup[track->getPoint(-1)];
if (kernel_from_event) if (kernel_from_event)
std::sort(available_obs.begin(), available_obs.end(), event_kernel_sort_fn); std::sort(available_obs.begin(), available_obs.end(), event_kernel_sort_fn);
...@@ -316,6 +321,8 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter ...@@ -316,6 +321,8 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter
std::sort(available_obs.begin(), available_obs.end(), backup_kernel_sort_fn); std::sort(available_obs.begin(), available_obs.end(), backup_kernel_sort_fn);
} }
// apply kernel value cutoff to avoid some difficult and
// unnecessary RK propagations
auto id = available_obs.back(); auto id = available_obs.back();
double kernel_val; double kernel_val;
if (kernel[current_id].count(id)) if (kernel[current_id].count(id))
...@@ -336,6 +343,7 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter ...@@ -336,6 +343,7 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter
track_point->setSortingParameter(seed_sort_param[id]); track_point->setSortingParameter(seed_sort_param[id]);
bool order_changed = track->sort(); bool order_changed = track->sort();
// fit hypotheses and remove bad hypotheses
std::vector<genfit::AbsTrackRep *> all_reps(track->getTrackReps()); std::vector<genfit::AbsTrackRep *> all_reps(track->getTrackReps());
std::vector<genfit::AbsTrackRep *> bad_reps; std::vector<genfit::AbsTrackRep *> bad_reps;
for (auto rep : all_reps) { for (auto rep : all_reps) {
...@@ -351,6 +359,9 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter ...@@ -351,6 +359,9 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter
if (pval > 0.001) continue; if (pval > 0.001) continue;
bad_reps.push_back(rep); bad_reps.push_back(rep);
} }
// if all hypotheses are bad, return the latest hit to the list of
// available hits and break
if (bad_reps.size() == all_reps.size()) { if (bad_reps.size() == all_reps.size()) {
for (int i = 1; i < track->getNumPoints(); i++) { for (int i = 1; i < track->getNumPoints(); i++) {
if (track_point == track->getPoint(i)) { if (track_point == track->getPoint(i)) {
...@@ -365,6 +376,8 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter ...@@ -365,6 +376,8 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter
for (auto rep : bad_reps) track->deleteTrackRep(track->getIdForRep(rep)); for (auto rep : bad_reps) track->deleteTrackRep(track->getIdForRep(rep));
} }
// determine best hypothesis
track->determineCardinalRep(); track->determineCardinalRep();
auto rep = track->getCardinalRep(); auto rep = track->getCardinalRep();
...@@ -474,17 +487,26 @@ int main(int argc, char **argv) { ...@@ -474,17 +487,26 @@ int main(int argc, char **argv) {
genfit::MaterialEffects::getInstance()->setNoEffects(); genfit::MaterialEffects::getInstance()->setNoEffects();
genfit::MaterialEffects::getInstance()->setEnergyLossBrems(false); genfit::MaterialEffects::getInstance()->setEnergyLossBrems(false);
genfit::KalmanFitter fitter; genfit::KalmanFitter fitter;
// set convergence conditions for track fitting
fitter.setDeltaPval(0.1); fitter.setDeltaPval(0.1);
fitter.setRelChi2Change(1); fitter.setRelChi2Change(1);
int comp = 0; int comp = 0;
int n_validate = 6; int n_seed = 4;
int n_validate = 8;
int opt; int opt;
while ((opt = getopt(argc, argv, "c:n:h")) != -1) { while ((opt = getopt(argc, argv, "c:s:n:h")) != -1) {
switch (opt) { switch (opt) {
case 'c': case 'c':
comp = atoi(optarg); comp = atoi(optarg);
break; break;
case 's':
n_seed = atoi(optarg);
if (n_seed < 1) {
std::cerr << "s must be >= 2" << std::endl;
exit(EXIT_FAILURE);
}
break;
case 'n': case 'n':
n_validate = atoi(optarg); n_validate = atoi(optarg);
if (n_validate < 2) { if (n_validate < 2) {
...@@ -521,12 +543,16 @@ int main(int argc, char **argv) { ...@@ -521,12 +543,16 @@ int main(int argc, char **argv) {
writer->SetCompression(proio::UNCOMPRESSED); writer->SetCompression(proio::UNCOMPRESSED);
} }
// loop over events
proio::Event event; proio::Event event;
int n_events = 0; int n_events = 0;
while (reader->Next(&event)) { while (reader->Next(&event)) {
std::cout << "event " << n_events << std::endl; std::cout << "event " << n_events << std::endl;
// if new field or geometry information is available, apply it
checkMetadata(&event); checkMetadata(&event);
trackEvent(&event, n_validate, &fitter); // fit tracks
trackEvent(&event, n_seed, n_validate, &fitter);
// write out modified event
writer->Push(&event); writer->Push(&event);
n_events++; n_events++;
} }
......
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