Added use of event kernel

src/EICTrack.cc

@@ -180,12 +180,34 @@ double forwardPVal(genfit::Track *track, genfit::AbsTrackRep *rep) {
 }
 
 void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter *fitter) {
+    // build kernel from event
+    std::map<uint64_t, std::map<uint64_t, double>> kernel;
+    bool kernel_from_event = false;
+    for (auto id : event->TaggedEntries("Reconstructed")) {
+        auto kernel_ud = dynamic_cast<eic::KernelMatrix *>(event->GetEntry(id));
+        if (kernel_ud) {
+            kernel_from_event = true;
+            int flat_index = 0;
+            for (int i = 0; i < kernel_ud->observation_size(); i++) {
+                uint64_t i_id = kernel_ud->observation(i);
+                auto &row = kernel[i_id];
+                for (int j = i; j < kernel_ud->observation_size(); j++) {
+                    uint64_t j_id = kernel_ud->observation(j);
+                    row[j_id] = kernel_ud->kflat(flat_index);
+                    flat_index++;
+                }
+            }
+            break;
+        }
+    }
+
     // build available observations list and a map of observations to
     // GenFit measurements
     std::vector<uint64_t> available_obs;
     std::map<uint64_t, genfit::AbsMeasurement *> measurements;
     std::map<uint64_t, double> seed_sort_param;
     for (auto id : event->TaggedEntries("Tracker")) {
+        if (kernel_from_event && kernel.count(id) == 0) continue;
         auto obs = dynamic_cast<eic::EnergyDep *>(event->GetEntry(id));
         if (obs) {
             auto measurement = makeMeasurement(obs);
@@ -200,9 +222,6 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter
                     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();
             }
         }
     }
@@ -216,18 +235,35 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter
         std::sort(available_obs.begin(), available_obs.end(), seed_sort_fn);
         track_obs.push_back(available_obs.back());
         available_obs.pop_back();
+        auto current_id = track_obs.back();
         auto measurement = measurements[track_obs[0]];
         TVector3 ip_pos(0, 0, 0);
         auto meas_coords = measurement->getRawHitCoords();
         TVector3 meas_pos(10 * meas_coords[0], 10 * meas_coords[1], 10 * meas_coords[2]);
-        auto 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;
+            if (kernel[current_id].count(a_id))
+                kernel_a = kernel[current_id][a_id];
+            else
+                kernel_a = kernel[a_id][current_id];
+            if (kernel[current_id].count(b_id))
+                kernel_b = kernel[current_id][b_id];
+            else
+                kernel_b = kernel[b_id][current_id];
+            return kernel_a < kernel_b;
+        };
+        auto backup_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]);
             auto b_coords = measurements[b_id]->getRawHitCoords();
             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();
         };
-        std::sort(available_obs.begin(), available_obs.end(), kernel_sort_fn);
+        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);
@@ -235,7 +271,7 @@ 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 = 3;
+        int n_seed = n_validate - 1;
         for (int i = 0; i < n_seed - 1; i++) {
             track_obs.push_back(available_obs.back());
             available_obs.pop_back();
@@ -261,7 +297,10 @@ void trackEvent(proio::Event *event, const int n_validate, genfit::KalmanFitter
         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]);
-            std::sort(available_obs.begin(), available_obs.end(), kernel_sort_fn);
+            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();