Write LOR

abracadabra/abracadabra.cc

@@ -5,6 +5,7 @@
 #include "geometries/nema.hh"
 #include "geometries/samples.hh"
 #include "geometries/sipm.hh"
+#include "utils/enumerate.hh"
 
 #include <G4RunManager.hh>
 #include <G4RunManagerFactory.hh>
@@ -25,6 +26,72 @@
 using std::make_unique;
 using std::unique_ptr;
 
+
+class write_with {
+public:
+  write_with(hdf5_io& writer)
+    : PROCESS_HITS{ [this](auto step){ return this -> process_hits(step); } }
+    , END_OF_EVENT{ [this](auto hc)  {        this -> end_of_event(hc  ); } }
+    , writer{writer}
+  {}
+
+  bool process_hits(G4Step* step) {
+    auto track    = step -> GetTrack();
+    auto particle = track -> GetParticleDefinition();
+    auto name     = particle -> GetParticleName();
+    if (name == "gamma") {
+      hits.push_back(*step);
+    }
+    return true;
+  }
+
+  void end_of_event(G4HCofThisEvent*) {
+    double event_id = n4::event_number();
+    double total_energy = 0;
+    std::vector<double> rs;
+    std::vector<double> phis;
+    std::vector<double> zs;
+    std::vector<double> ts;
+
+    // Save only events where two gammas has been detected
+    if (hits.size() == 2){
+      std::cout << "Two gammas detected!" << std::endl;
+
+      for (auto [i, step]: enumerate(hits)){
+          auto pos      = step . GetPostStepPoint() -> GetPosition();
+          auto track    = step . GetTrack();
+          auto energy   = track -> GetKineticEnergy();
+          auto particle = track -> GetParticleDefinition();
+          auto name     = particle -> GetParticleName();
+          auto id       = track -> GetTrackID();
+          auto time     = track -> GetGlobalTime();
+
+          auto r   = sqrt( pos.x()*pos.x() + pos.y()*pos.y() );
+          auto phi = atan2( pos.y(), pos.x() );
+
+          total_energy += energy;
+          rs  .push_back(r);
+          phis.push_back(phi);
+          zs  .push_back(pos.z());
+          ts  .push_back(time);
+      }
+      writer.write_lor_info(event_id, total_energy, rs[0], phis[0], zs[0], ts[0], rs[1], phis[1], zs[1], ts[1]);
+    }
+    hits = {};
+  }
+
+  // TODO: better docs: use for filling slots in n4::sensitive_detector
+  n4::sensitive_detector::process_hits_fn const PROCESS_HITS;
+  n4::sensitive_detector::end_of_event_fn const END_OF_EVENT;
+
+private:
+  std::vector<G4Step> hits{};
+  hdf5_io& writer;
+
+};
+
+
+
 int main(int argc, char** argv) {
   // Detect interactive mode (if no arguments) and define UI session
   auto ui = argc == 1
@@ -39,7 +106,18 @@ int main(int argc, char** argv) {
     {G4RunManagerFactory::CreateRunManager(G4RunManagerType::Serial)};
 
   // For use with phantom_in_cylinder
-  auto phantom = a_nema_phantom();
+  auto phantom = build_nema_phantom{}
+    .activity(0)
+    .length(140*mm)
+    .inner_radius(58*mm)
+    .outer_radius(76*mm)
+    .sphere( 5  *mm / 2, 20)
+    .sphere( 6.5*mm / 2, 20)
+    .sphere( 8.5*mm / 2, 20)
+    .sphere(11  *mm / 2, 20)
+    .sphere(14  *mm / 2, 0)
+    .sphere(19  *mm / 2, 0)
+    .build();
 
   int last_interesting_event = 0;
   size_t count_interesting_event = 0;
@@ -105,8 +183,15 @@ int main(int argc, char** argv) {
     got_electron = false;
   };
 
+
+  // Add writer for LORs
+  // TODO: Filename should be taken from config file
+  std::string hdf5_file_name = "test_lor.h5";
+  auto writer = new hdf5_io{hdf5_file_name};
+  auto fwd = write_with{*writer};
+
   // pick one:
-  auto sd = new n4::sensitive_detector{"Noisy_detector", make_noise, eoe};
+  auto sd = new n4::sensitive_detector{"Noisy_detector", fwd.PROCESS_HITS, fwd.END_OF_EVENT};
   //auto sd = nullptr;
 
   // Set mandatory initialization classes
@@ -177,11 +262,11 @@ int main(int argc, char** argv) {
         ui_manager->ApplyCommand("/vis/viewer/set/viewpointThetaPhi "
                                  + std::to_string(theta) + ' ' + std::to_string(phi));
       };
-      {
-        nain4::silence _{G4cout};
-        for (int phi  =PHI  ; phi  <360+PHI   ; phi  +=4) { view(THETA, phi); }
-        for (int theta=THETA; theta<360+THETAF; theta+=4) { view(theta, PHI); }
-      }
+      // {
+      //   nain4::silence _{G4cout};
+      //   for (int phi  =PHI  ; phi  <360+PHI   ; phi  +=4) { view(THETA, phi); }
+      //   for (int theta=THETA; theta<360+THETAF; theta+=4) { view(theta, PHI); }
+      // }
       ui -> SessionStart();
     }
   }

abracadabra/src/geometries/samples.cc

@@ -56,7 +56,7 @@ G4VPhysicalVolume* cylinder_lined_with_hamamatsus(double length, double radius,
 
   auto cavity_r = radius - dr_Xe;
 
-  auto xenon    = volume<G4Tubs>("LXe"     , lXe, cavity_r,   radius, length/2, 0.0, CLHEP::twopi);
+  auto xenon    = volume<G4Tubs>("LXe"     , air, cavity_r,   radius, length/2, 0.0, CLHEP::twopi);
   auto cavity   = volume<G4Tubs>("Cavity"  , air, 0.0     , cavity_r, length/2, 0.0, CLHEP::twopi);
   auto envelope = volume<G4Box> ("Envelope", air, 1.1*radius, 1.1*radius, 1.1*length/2);
 

abracadabra/src/geometries/samples.hh

@@ -5,6 +5,6 @@
 
 nema_phantom       a_nema_phantom();
 G4VPhysicalVolume* square_array_of_sipms(G4VSensitiveDetector* = nullptr);
-G4VPhysicalVolume* phantom_in_cylinder(nema_phantom const&, double, double, G4VSensitiveDetector* = nullptr);
+G4VPhysicalVolume* phantom_in_cylinder(nema_phantom const&, G4double length, G4double dr_Xe, G4VSensitiveDetector* = nullptr);
 G4VPhysicalVolume* cylinder_lined_with_hamamatsus(double length, double radius, double dr_Xe,
                                                   G4VSensitiveDetector* = nullptr);

abracadabra/src/io/hdf5.cc

@@ -9,7 +9,8 @@ namespace HF { using namespace HighFive; }
 hdf5_io::hdf5_io(std::string fname)
 : filename{fname}
 , runinfo_index{0}
-, hit_index{0} {}
+, hit_index{0}
+, lor_index{0} {}
 
 
 HF::CompoundType create_hit_type() {
@@ -27,6 +28,33 @@ HF::CompoundType create_runinfo_type() {
 }
 HIGHFIVE_REGISTER_TYPE(run_info_t, create_runinfo_type)
 
+
+HF::CompoundType create_lor_type() {
+  return {{"event_id"   , HF::AtomicType<double>{}},
+          {"true_energy", HF::AtomicType<double>{}},
+          {"true_r1"    , HF::AtomicType<double>{}},
+          {"true_phi1"  , HF::AtomicType<double>{}},
+          {"true_z1"    , HF::AtomicType<double>{}},
+          {"true_t1"    , HF::AtomicType<double>{}},
+          {"true_r2"    , HF::AtomicType<double>{}},
+          {"true_phi2"  , HF::AtomicType<double>{}},
+          {"true_z2"    , HF::AtomicType<double>{}},
+          {"true_t2"    , HF::AtomicType<double>{}},
+          {"phot_like1" , HF::AtomicType<double>{}},
+          {"phot_like2" , HF::AtomicType<double>{}},
+          {"reco_r1"    , HF::AtomicType<double>{}},
+          {"reco_phi1"  , HF::AtomicType<double>{}},
+          {"reco_z1"    , HF::AtomicType<double>{}},
+          {"reco_t1"    , HF::AtomicType<double>{}},
+          {"reco_r2"    , HF::AtomicType<double>{}},
+          {"reco_phi2"  , HF::AtomicType<double>{}},
+          {"reco_z2"    , HF::AtomicType<double>{}},
+          {"reco_t2"    , HF::AtomicType<double>{}},
+          {"not_sel"    , HF::AtomicType<double>{}}};
+}
+HIGHFIVE_REGISTER_TYPE(lor_t, create_lor_type)
+
+
 void set_string_param(char * to, const char * from, unsigned int max_len) {
   memset(to, 0, max_len);
   strcpy(to, from);
@@ -53,6 +81,9 @@ void hdf5_io::ensure_open_for_writing() {
 
   group.createDataSet("hits"         , dataspace, create_hit_type()    , props);
   group.createDataSet("configuration", dataspace, create_runinfo_type(), props);
+
+  group = open_for_writing->createGroup("reco_info");
+  group.createDataSet("table", dataspace, create_lor_type()    , props);
 }
 
 void hdf5_io::write_run_info(const char* param_key, const char* param_value) {
@@ -95,3 +126,44 @@ std::vector<hit_t> hdf5_io::read_hit_info() {
   hits_table.read(hits);
   return hits;
 }
+
+
+void hdf5_io::write_lor_info(double event_id, double energy,
+          double r1, double phi1, double z1, double t1,
+          double r2, double phi2, double z2, double t2){
+  lor_t lor;
+  lor.event_id    = event_id;
+  lor.true_energy = energy;
+  lor.true_r1     = r1;
+  lor.true_phi1   = phi1;
+  lor.true_z1     = z1;
+  lor.true_t1     = t1;
+  lor.true_r2     = r2;
+  lor.true_phi2   = phi2;
+  lor.true_z2     = z2;
+  lor.true_t2     = t2;
+  lor.phot_like1  = 1;
+  lor.phot_like2  = 1;
+  lor.reco_r1     = r1;
+  lor.reco_phi1   = phi1;
+  lor.reco_z1     = z1;
+  lor.reco_t1     = t1;
+  lor.reco_r2     = r2;
+  lor.reco_phi2   = phi2;
+  lor.reco_z2     = z2;
+  lor.reco_t2     = t2;
+  lor.not_sel     = 0;
+
+  std::vector<lor_t> data{lor};
+  unsigned int n_elements = data.size();
+
+  ensure_open_for_writing();
+  HF::Group   group = open_for_writing -> getGroup("reco_info");
+  HF::DataSet table = group.getDataSet("table");
+
+  // Create extra space in the table and append the new data
+  table.resize({lor_index + n_elements});
+  table.select({lor_index}, {n_elements}).write(data);
+
+  lor_index += n_elements;
+}

abracadabra/src/io/hdf5.hh

@@ -19,6 +19,30 @@ typedef struct {
   double t;
 } hit_t;
 
+typedef struct {
+  double event_id;
+  double true_energy;
+  double true_r1;
+  double true_phi1;
+  double true_z1;
+  double true_t1;
+  double true_r2;
+  double true_phi2;
+  double true_z2;
+  double true_t2;
+  double phot_like1;
+  double phot_like2;
+  double reco_r1;
+  double reco_phi1;
+  double reco_z1;
+  double reco_t1;
+  double reco_r2;
+  double reco_phi2;
+  double reco_z2;
+  double reco_t2;
+  double not_sel;
+} lor_t;
+
 
 // TODO There's something fishy about the implementation behind this interface:
 // it holds on to a filename, and then opens the file each time it wants to
@@ -30,6 +54,9 @@ public:
 
   void write_run_info(const char* param_key, const char* param_value);
   void write_hit_info(unsigned int evt_id, double x, double y, double z, double t);
+  void write_lor_info(double event_id, double energy,
+          double r1, double phi1, double z1, double t1,
+          double r2, double phi2, double z2, double t2);
   void flush() { if (open_for_writing) { open_for_writing -> flush(); } }
 
   std::vector<hit_t> read_hit_info();
@@ -42,6 +69,7 @@ private:
   std::string filename;
   unsigned int runinfo_index;
   unsigned int hit_index;
+  unsigned int lor_index;
   std::optional<HighFive::File> open_for_writing;
 };
 
@@ -51,5 +79,4 @@ typedef struct {
   char param_value[hdf5_io::CONFLEN];
 } run_info_t;
 
-
 #endif