diff --git a/abracadabra/abracadabra.cc b/abracadabra/abracadabra.cc index 0a96f8e..8416b39 100644 --- a/abracadabra/abracadabra.cc +++ b/abracadabra/abracadabra.cc @@ -5,6 +5,7 @@ #include "geometries/nema.hh" #include "geometries/samples.hh" #include "geometries/sipm.hh" +#include "utils/enumerate.hh" #include #include @@ -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 rs; + std::vector phis; + std::vector zs; + std::vector 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 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(); } } diff --git a/abracadabra/src/geometries/samples.cc b/abracadabra/src/geometries/samples.cc index dfd0a4a..7125ffc 100644 --- a/abracadabra/src/geometries/samples.cc +++ b/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("LXe" , lXe, cavity_r, radius, length/2, 0.0, CLHEP::twopi); + auto xenon = volume("LXe" , air, cavity_r, radius, length/2, 0.0, CLHEP::twopi); auto cavity = volume("Cavity" , air, 0.0 , cavity_r, length/2, 0.0, CLHEP::twopi); auto envelope = volume ("Envelope", air, 1.1*radius, 1.1*radius, 1.1*length/2); diff --git a/abracadabra/src/geometries/samples.hh b/abracadabra/src/geometries/samples.hh index e2bbc44..21ae326 100644 --- a/abracadabra/src/geometries/samples.hh +++ b/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); diff --git a/abracadabra/src/io/hdf5.cc b/abracadabra/src/io/hdf5.cc index de38d08..6e92e7a 100644 --- a/abracadabra/src/io/hdf5.cc +++ b/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{}}, + {"true_energy", HF::AtomicType{}}, + {"true_r1" , HF::AtomicType{}}, + {"true_phi1" , HF::AtomicType{}}, + {"true_z1" , HF::AtomicType{}}, + {"true_t1" , HF::AtomicType{}}, + {"true_r2" , HF::AtomicType{}}, + {"true_phi2" , HF::AtomicType{}}, + {"true_z2" , HF::AtomicType{}}, + {"true_t2" , HF::AtomicType{}}, + {"phot_like1" , HF::AtomicType{}}, + {"phot_like2" , HF::AtomicType{}}, + {"reco_r1" , HF::AtomicType{}}, + {"reco_phi1" , HF::AtomicType{}}, + {"reco_z1" , HF::AtomicType{}}, + {"reco_t1" , HF::AtomicType{}}, + {"reco_r2" , HF::AtomicType{}}, + {"reco_phi2" , HF::AtomicType{}}, + {"reco_z2" , HF::AtomicType{}}, + {"reco_t2" , HF::AtomicType{}}, + {"not_sel" , HF::AtomicType{}}}; +} +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 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 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; +} diff --git a/abracadabra/src/io/hdf5.hh b/abracadabra/src/io/hdf5.hh index 91505ec..968740e 100644 --- a/abracadabra/src/io/hdf5.hh +++ b/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 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 open_for_writing; }; @@ -51,5 +79,4 @@ typedef struct { char param_value[hdf5_io::CONFLEN]; } run_info_t; - #endif diff --git a/abracadabra/test_lor.h5 b/abracadabra/test_lor.h5 new file mode 100644 index 0000000..3e72a18 Binary files /dev/null and b/abracadabra/test_lor.h5 differ