Merge lp:~jan.greis/maus/global-recon-rk into lp:maus/merge

[Jan Greis]

=== modified file 'src/common_cpp/Recon/Global/ImportKLRecon.cc'
--- src/common_cpp/Recon/Global/ImportKLRecon.cc	2016-08-31 11:59:03 +0000
+++ src/common_cpp/Recon/Global/ImportKLRecon.cc	2016-10-19 13:41:35 +0000
@@ -25,33 +25,84 @@
 namespace recon {
 namespace global {
 
-  void ImportKLRecon::process(const MAUS::KLEvent &kl_event,
-			       MAUS::GlobalEvent* global_event,
-			       std::string mapper_name) {
-
-    KLEventCellHit KLEventCellHit = kl_event.GetKLEventCellHit();
-
-    for (size_t i = 0; i < KLEventCellHit.GetKLCellHitArraySize(); ++i) {
-      KLCellHit KLCellHit = KLEventCellHit.GetKLCellHitArrayElement(i);
-      MAUS::DataStructure::Global::SpacePoint* GlobalKLSpacePoint =
-	new MAUS::DataStructure::Global::SpacePoint();
-      GlobalKLSpacePoint->set_ADC_charge(KLCellHit.GetCharge());
-      GlobalKLSpacePoint->set_ADC_charge_product(KLCellHit.GetChargeProduct());
-      GlobalKLSpacePoint->set_detector(MAUS::DataStructure::Global::kCalorimeter);
-      double x = KLCellHit.GetGlobalPosX();
-      double x_err = KLCellHit.GetErrorX();
-      double y = KLCellHit.GetGlobalPosY();
-      double y_err = KLCellHit.GetErrorY();
-      double z = KLCellHit.GetGlobalPosZ();
-      double z_err = KLCellHit.GetErrorZ();
-      double t = -1000000;
-      double t_err = 1000000;
+  void ImportKLRecon::process(const MAUS::KLEvent &kl_event, MAUS::GlobalEvent* global_event,
+      std::string mapper_name, bool merge_cell_hits) {
+
+    KLEventCellHit kl_event_cell_hit = kl_event.GetKLEventCellHit();
+    std::vector<KLCellHit> kl_cell_hits = kl_event_cell_hit.GetKLCellHitArray();
+    if (kl_cell_hits.size() < 1) {
+      return;
+    }
+    std::sort(kl_cell_hits.begin(), kl_cell_hits.end(), KLCellHitSort);
+    std::vector<std::vector<size_t> > cell_groupings = AdjacentCellGroupings(kl_cell_hits,
+        merge_cell_hits);
+
+    for (size_t i = 0; i < cell_groupings.size(); i++) {
+      double chargesum = 0.0, chargeproductsum = 0.0;
+      double x = 0.0, x_err = 0.0, y = 0.0, y_err = 0.0, z = 0.0, z_err = 0.0;
+      double weighted_y2 = 0.0;
+      double t = -1000000, t_err = 1000000;
+      for (size_t j = 0; j < cell_groupings.at(i).size(); j++) {
+        KLCellHit kl_cell_hit = kl_cell_hits.at(cell_groupings.at(i).at(j));
+        double charge = kl_cell_hit.GetCharge();
+        double chargeproduct = kl_cell_hit.GetChargeProduct();
+        chargesum += charge;
+        chargeproductsum += chargeproduct;
+        x += kl_cell_hit.GetGlobalPosX()*charge;
+        y += kl_cell_hit.GetGlobalPosY()*charge;
+        z += kl_cell_hit.GetGlobalPosZ()*charge;
+        weighted_y2 += kl_cell_hit.GetGlobalPosY()*kl_cell_hit.GetGlobalPosY()*charge;
+        // Are the same for all, so this way we can easily just pick out the last one
+        x_err = kl_cell_hit.GetErrorX();
+        z_err = kl_cell_hit.GetErrorZ();
+        if (cell_groupings.at(i).size() == 1) {
+          y_err = kl_cell_hit.GetErrorY();
+        }
+      }
+      x /= chargesum;
+      y /= chargesum;
+      z /= chargesum;
+      if (cell_groupings.at(i).size() != 1) {
+        weighted_y2 /= chargesum;
+        y_err = std::sqrt(weighted_y2 - y*y);
+      }
       TLorentzVector pos(x, y, z, t);
       TLorentzVector pos_err(x_err, y_err, z_err, t_err);
-      GlobalKLSpacePoint->set_position(pos);
-      GlobalKLSpacePoint->set_position_error(pos_err);
-      global_event->add_space_point(GlobalKLSpacePoint);
-    }
+      DataStructure::Global::SpacePoint* global_space_point =
+            new DataStructure::Global::SpacePoint();
+      global_space_point->set_ADC_charge(chargesum);
+      global_space_point->set_ADC_charge_product(chargeproductsum);
+      global_space_point->set_detector(MAUS::DataStructure::Global::kCalorimeter);
+      global_space_point->set_position(pos);
+      global_space_point->set_position_error(pos_err);
+      global_event->add_space_point(global_space_point);
+    }
+  }
+
+  std::vector<std::vector<size_t> > ImportKLRecon::AdjacentCellGroupings(
+      std::vector<KLCellHit> kl_cell_hits, bool merge_cell_hits) {
+    std::vector<std::vector<size_t> > cell_groupings;
+    std::vector<size_t> cell_hit_ids;
+    if (kl_cell_hits.size() == 1) {
+      cell_hit_ids.push_back(0);
+      cell_groupings.push_back(cell_hit_ids);
+    } else {
+      for (size_t i = 0; i < kl_cell_hits.size() - 1; i++) {
+        cell_hit_ids.push_back(i);
+        if (kl_cell_hits.at(i+1).GetCell() != kl_cell_hits.at(i).GetCell() + 1
+            or !merge_cell_hits) { // If configuration is set not to merge, one std::vector for each
+          cell_groupings.push_back(cell_hit_ids);
+          cell_hit_ids.clear();
+        }
+      }
+      cell_hit_ids.push_back(kl_cell_hits.size() - 1);
+      cell_groupings.push_back(cell_hit_ids);
+    }
+    return cell_groupings;
+  }
+
+  bool ImportKLRecon::KLCellHitSort(KLCellHit hit1, KLCellHit hit2) {
+    return (hit1.GetCell() < hit2.GetCell());
   }
 
 } // ~namespace global

=== modified file 'src/common_cpp/Recon/Global/ImportKLRecon.hh'
--- src/common_cpp/Recon/Global/ImportKLRecon.hh	2014-06-20 10:51:38 +0000
+++ src/common_cpp/Recon/Global/ImportKLRecon.hh	2016-10-19 13:41:35 +0000
@@ -47,14 +47,26 @@
     /// Destructor
     ~ImportKLRecon() {}
 
-    /**  Main process, accepting the MAUS::KLEvent and importing
+    /** @brief Main process, accepting the MAUS::KLEvent and importing
      *   space points into an existing MAUS::GlobalEvent
      *  @param kl_event  The reconstructed KL Event
      *  @param global_event The Global Event, which will be changed
      */
-    void process(const MAUS::KLEvent &kl_event,
-		 MAUS::GlobalEvent* global_event,
-		 std::string mapper_name);
+    void process(const MAUS::KLEvent &kl_event, MAUS::GlobalEvent* global_event,
+                 std::string mapper_name, bool merge_cell_hits);
+
+    /** @brief Determine whether groupings of adjacent cells exist in the KLCellHits.
+     *    Returns a vector of vectors of cell numbers appropriately grouped, e.g. [[4],[7,8]]
+     *    If merge_cell_hits is false, all are separated, e.g. [[4],[7],[8]]
+     */
+    std::vector<std::vector<size_t> > AdjacentCellGroupings(std::vector<KLCellHit> kl_cell_hits,
+        bool merge_cell_hits);
+
+    /**
+     * @brief Returns whether the first KLCellHit has a lower cell number than the
+     * second one. Used to perform an std::sort on a vector of KLCellHits.
+     */
+    static bool KLCellHitSort(KLCellHit hit1, KLCellHit hit2);
 
   private:
     /// Disallow copy constructor as unnecessary

=== modified file 'src/common_py/ConfigurationDefaults.py'
--- src/common_py/ConfigurationDefaults.py	2016-09-13 09:30:18 +0000
+++ src/common_py/ConfigurationDefaults.py	2016-10-19 13:41:35 +0000
@@ -826,3 +826,7 @@
   "TOF2":0.0,
   "KL":0.0
 }
+
+# Whether multiple adjacent cell hits in the KL should be merged into single spacepoints on import
+# into the global datastructure
+global_merge_kl_cell_hits = True

=== modified file 'src/map/MapCppGlobalReconImport/MapCppGlobalReconImport.cc'
--- src/map/MapCppGlobalReconImport/MapCppGlobalReconImport.cc	2016-08-31 11:59:03 +0000
+++ src/map/MapCppGlobalReconImport/MapCppGlobalReconImport.cc	2016-10-19 13:41:35 +0000
@@ -55,6 +55,7 @@
 			    "MapCppGlobalReconImport::_birth");
     }
 
+    _merge_cell_hits = _configJSON["global_merge_kl_cell_hits"].asBool();
     _configCheck = true;
     _classname = "MapCppGlobalReconImport";
   }
@@ -113,7 +114,7 @@
 	if (recon_event->GetKLEvent()) {
 	  MAUS::KLEvent* kl_event = recon_event->GetKLEvent();
 	  MAUS::recon::global::ImportKLRecon klrecon_importer;
-	  klrecon_importer.process((*kl_event), global_event, _classname);
+	  klrecon_importer.process((*kl_event), global_event, _classname, _merge_cell_hits);
 	}
 	// Import EMR tracks into global event
 	if (recon_event->GetEMREvent()) {

=== modified file 'src/map/MapCppGlobalReconImport/MapCppGlobalReconImport.hh'
--- src/map/MapCppGlobalReconImport/MapCppGlobalReconImport.hh	2016-08-31 11:59:03 +0000
+++ src/map/MapCppGlobalReconImport/MapCppGlobalReconImport.hh	2016-10-19 13:41:35 +0000
@@ -95,6 +95,8 @@
     Json::Value _configJSON;
     ///  JsonCpp setup
     Json::Reader _reader;
+    /// Whether adjacent cell hits should be merged into single spacepoints
+    bool _merge_cell_hits;
     /// Mapper name, useful for tracking results...
     std::string _classname;
     // Geometry file name