From ee1c0e69847bece5a0416362798ba1ebcf129128 Mon Sep 17 00:00:00 2001 From: Stephen Fegan Date: Thu, 4 Jul 2024 16:35:56 +0200 Subject: [PATCH] Modify GroundMap to handle multiple detectors --- include/simulation/vcl_iact_ground_map.hpp | 175 ++++++++++++------- proto/simulation/vcl_iact.proto | 38 ++++ unit_tests/simulation/speedtest_vcl_iact.cpp | 4 +- 3 files changed, 147 insertions(+), 70 deletions(-) diff --git a/include/simulation/vcl_iact_ground_map.hpp b/include/simulation/vcl_iact_ground_map.hpp index e88c3a3e..b9dda2a0 100644 --- a/include/simulation/vcl_iact_ground_map.hpp +++ b/include/simulation/vcl_iact_ground_map.hpp @@ -26,7 +26,9 @@ #include +#include #include +#include namespace calin { namespace simulation { namespace vcl_iact { @@ -52,31 +54,35 @@ template class alignas(VCLArchitecture::vec_bytes) VCL using double_at = typename VCLArchitecture::double_at; #endif // not defined SWIG - VCLIACTGroundMap(calin::simulation::atmosphere::LayeredRefractiveAtmosphere* atm, double dzatm_profile, - const calin::ix::simulation::vcl_iact::VCLIACTConfiguration& config = VCLIACTTrackVisitor::default_config(), - calin::math::rng::VCLRNG* rng = nullptr, - bool adopt_atm = false, bool adopt_rng = false); VCLIACTGroundMap(calin::simulation::atmosphere::LayeredRefractiveAtmosphere* atm, - const calin::ix::simulation::vcl_iact::VCLIACTConfiguration& config = VCLIACTTrackVisitor::default_config(), + const calin::ix::simulation::vcl_iact::VCLIACTGroundMapConfiguration& config = default_config(), calin::math::rng::VCLRNG* rng = nullptr, bool adopt_atm = false, bool adopt_rng = false); virtual ~VCLIACTGroundMap(); + static calin::ix::simulation::vcl_iact::VCLIACTGroundMapConfiguration default_config(); + + const calin::ix::simulation::vcl_iact::VCLIACTGroundMapConfiguration& config() const { return config_; } double dzatm_profile() const { return 1.0/dzatm_profile_inv_; } - double ground_radius_cut() const { return std::sqrt(r2gnd_cut_); } - void set_ground_radius_cut(double r) { r2gnd_cut_ = r*r; } - void set_store_emission_point(bool store = true) { store_emission_pt_ = store; } - - const std::vector& xatm() const { return xatm_; } - const std::vector& yatm() const { return yatm_; } - const std::vector& zatm() const { return zatm_; } - const std::vector& tgnd() const { return tgnd_; } - const std::vector& xgnd() const { return xgnd_; } - const std::vector& ygnd() const { return ygnd_; } - const std::vector& uxgnd() const { return uxgnd_; } - const std::vector& uygnd() const { return uygnd_; } - double ncherenkov() const { return ncherenkov_; } + const std::vector& xatm(unsigned idetector) const { + return detector_.at(idetector).xatm; } + const std::vector& yatm(unsigned idetector) const { + return detector_.at(idetector).yatm; } + const std::vector& zatm(unsigned idetector) const { + return detector_.at(idetector).zatm; } + const std::vector& tgnd(unsigned idetector) const { + return detector_.at(idetector).tgnd; } + const std::vector& xgnd(unsigned idetector) const { + return detector_.at(idetector).xgnd; } + const std::vector& ygnd(unsigned idetector) const { + return detector_.at(idetector).ygnd; } + const std::vector& uxgnd(unsigned idetector) const { + return detector_.at(idetector).uxgnd; } + const std::vector& uygnd(unsigned idetector) const { + return detector_.at(idetector).uygnd; } + + double num_cherenkov() const { return ncherenkov_; } const std::vector& zatm_profile() const { return zatm_profile_; } #ifndef SWIG @@ -85,23 +91,39 @@ template class alignas(VCLArchitecture::vec_bytes) VCL double_vt bandwidth, double_vt ray_weight) final; protected: + struct Detector { + double r2; + const calin::ix::simulation::vcl_iact::VCLIACTGroundMapDetectorConfiguration* config; + std::vector xatm; + std::vector yatm; + std::vector zatm; + std::vector tgnd; + std::vector xgnd; + std::vector ygnd; + std::vector uxgnd; + std::vector uygnd; + void clear() { + xatm.clear(); + yatm.clear(); + zatm.clear(); + tgnd.clear(); + xgnd.clear(); + ygnd.clear(); + uxgnd.clear(); + uygnd.clear(); + } + }; + + calin::ix::simulation::vcl_iact::VCLIACTGroundMapConfiguration config_; + double ncherenkov_; std::vector zatm_profile_; double dzatm_profile_inv_; - std::vector xatm_; - std::vector yatm_; - std::vector zatm_; - std::vector tgnd_; - std::vector xgnd_; - std::vector ygnd_; - std::vector uxgnd_; - std::vector uygnd_; + std::vector detector_; unsigned iobs_ = 0; double zobs_ = 0.0; - double r2gnd_cut_ = 0.0; - bool store_emission_pt_ = false; #endif }; @@ -110,33 +132,45 @@ template class alignas(VCLArchitecture::vec_bytes) VCL template VCLIACTGroundMap:: VCLIACTGroundMap( calin::simulation::atmosphere::LayeredRefractiveAtmosphere* atm, - double dzatm_profile, - const calin::ix::simulation::vcl_iact::VCLIACTConfiguration& config, + const calin::ix::simulation::vcl_iact::VCLIACTGroundMapConfiguration& config, calin::math::rng::VCLRNG* rng, bool adopt_atm, bool adopt_rng): - VCLIACTTrackVisitor(atm, config, rng, adopt_atm, adopt_rng), - dzatm_profile_inv_(1.0/dzatm_profile) + VCLIACTTrackVisitor(atm, config.base_config(), rng, adopt_atm, adopt_rng), + config_(config), dzatm_profile_inv_(1.0/config.dzatm_profile()), + iobs_(config.observation_level()), zobs_(atm->zobs(iobs_)) { - iobs_ = 0; - zobs_ = this->atm_->zobs(iobs_); + using calin::math::special::SQR; + zatm_profile_.resize(std::ceil(atm->top_of_atmosphere() * dzatm_profile_inv_)); + detector_.resize(config_.detector_size()); + for(int idetector=0;idetector VCLIACTGroundMap:: -VCLIACTGroundMap( - calin::simulation::atmosphere::LayeredRefractiveAtmosphere* atm, - const calin::ix::simulation::vcl_iact::VCLIACTConfiguration& config, - calin::math::rng::VCLRNG* rng, - bool adopt_atm, bool adopt_rng): - VCLIACTGroundMap(atm, 5000.0 /* 50m */, config, rng, adopt_atm, adopt_rng) +~VCLIACTGroundMap() { // nothing to see here } -template VCLIACTGroundMap:: -~VCLIACTGroundMap() +template +calin::ix::simulation::vcl_iact::VCLIACTGroundMapConfiguration +VCLIACTGroundMap::default_config() { - // nothing to see here + calin::ix::simulation::vcl_iact::VCLIACTGroundMapConfiguration config; + config.mutable_base_config()->CopyFrom(VCLIACTTrackVisitor::default_config()); + auto* detector = config.add_detector(); + detector->set_x_gnd(0); + detector->set_y_gnd(0); + detector->set_r_gnd(1e5 /* 1km */); + detector->set_store_position(true); + detector->set_store_direction(true); + detector->set_store_time(true); + detector->set_store_emission_point(false); + detector->set_store_fraction(1.0); + return config; } template void VCLIACTGroundMap:: @@ -144,15 +178,9 @@ visit_event(const calin::simulation::tracker::Event& event, bool& kill_event) { ncherenkov_ = 0.0; std::fill(zatm_profile_.begin(), zatm_profile_.end(), 0.0); - xatm_.clear(); - yatm_.clear(); - zatm_.clear(); - tgnd_.clear(); - xgnd_.clear(); - ygnd_.clear(); - uxgnd_.clear(); - uygnd_.clear(); - + for(auto& idetector : detector_) { + idetector.clear(); + } return VCLIACTTrackVisitor::visit_event(event, kill_event); } @@ -185,10 +213,6 @@ propagate_rays(calin::math::ray::VCLRay ray, double_bvt ray_mask, ray_mask = ray.propagate_to_z_plane_with_mask(ray_mask, zobs_, false); - if(r2gnd_cut_) { - ray_mask &= ray.x()*ray.x() + ray.y()*ray.y() < r2gnd_cut_; - } - double_at t; double_at x; double_at y; @@ -200,19 +224,34 @@ propagate_rays(calin::math::ray::VCLRay ray, double_bvt ray_mask, ray.ux().store(ux); ray.uy().store(uy); - for(unsigned iv=0; iv::rng_->uniform_double(); + + for(unsigned idetector=0;idetectorx_gnd(); + double_vt yrel = ray.y() - detector.config->y_gnd(); + double_bvt store_mask = ray_mask && (xrel*xrel + yrel*yrel)store_fraction(); + + for(unsigned iv=0; ivstore_position()) { + detector.xgnd.push_back(x[iv]); + detector.ygnd.push_back(y[iv]); + } + if(detector.config->store_direction()) { + detector.uxgnd.push_back(ux[iv]); + detector.uygnd.push_back(uy[iv]); + } + if(detector.config->store_time()) { + detector.tgnd.push_back(t[iv]); + } + if(detector.config->store_time()) { + detector.xatm.push_back(atm_x[iv]); + detector.yatm.push_back(atm_y[iv]); + detector.zatm.push_back(atm_z[iv]); + } + } } - tgnd_.push_back(t[iv]); - xgnd_.push_back(x[iv]); - ygnd_.push_back(y[iv]); - uxgnd_.push_back(ux[iv]); - uygnd_.push_back(uy[iv]); - } } } diff --git a/proto/simulation/vcl_iact.proto b/proto/simulation/vcl_iact.proto index 0def7db3..b6341a05 100644 --- a/proto/simulation/vcl_iact.proto +++ b/proto/simulation/vcl_iact.proto @@ -84,3 +84,41 @@ message VCLIACTArrayConfiguration RefractionMode refraction_mode = 30 [(CFO).desc = "Atmospheric refraction mode." ]; }; + +message VCLIACTGroundMapDetectorConfiguration +{ + double x_gnd = 1 + [(CFO).desc = "X-osition of center of detector on ground.", + (CFO).units = "cm" ]; + double y_gnd = 2 + [(CFO).desc = "Y-position of center of detector on ground.", + (CFO).units = "cm" ]; + double r_gnd = 3 + [(CFO).desc = "Radius of detector on ground.", + (CFO).units = "cm" ]; + bool store_position = 4 + [(CFO).desc = "Store the position of photons in this detector." ]; + bool store_direction = 5 + [(CFO).desc = "Store the dirction of photons in this detector." ]; + bool store_time = 6 + [(CFO).desc = "Store the dirction of photons in this detector." ]; + bool store_emission_point = 7 + [(CFO).desc = "Store the emission point of the photons in this detector." ]; + double store_fraction = 8 + [(CFO).desc = "Fraction of photons to store." ]; +}; + +message VCLIACTGroundMapConfiguration +{ + VCLIACTConfiguration base_config = 1 + [(CFO).desc = "Base configuration." ]; + repeated VCLIACTGroundMapDetectorConfiguration detector = 2 + [(CFO).desc = "Ground detectors." ]; + double dzatm_profile = 3 + [(CFO).desc = "Width of bins in atmospheric profile.", + (CFO).units = "cm" ]; + uint32 observation_level = 4 + [(CFO).desc = "In the current implementation all telescopes must belong to " + "a single observation level." ]; +}; + diff --git a/unit_tests/simulation/speedtest_vcl_iact.cpp b/unit_tests/simulation/speedtest_vcl_iact.cpp index a2feaac9..c6b9b6d6 100644 --- a/unit_tests/simulation/speedtest_vcl_iact.cpp +++ b/unit_tests/simulation/speedtest_vcl_iact.cpp @@ -44,8 +44,8 @@ TEST(SpeedTestVCLIACT, Generate100Protons1TeV) { auto config = calin::simulation::vcl_iact::VCLIACTGroundMap< calin::util::vcl::VCL256Architecture>::default_config(); if(global_argc > 1) { - config.set_bandwidth(calin::util::string::double_from_string(global_argv[1])); - std::cerr << "HELLO: " << global_argv[1] << ' ' << config.bandwidth() << '\n'; + config.mutable_base_config()->set_bandwidth(calin::util::string::double_from_string(global_argv[1])); + std::cerr << "HELLO: " << global_argv[1] << ' ' << config.base_config().bandwidth() << '\n'; } auto* act = new calin::simulation::vcl_iact::VCLIACTGroundMap< calin::util::vcl::VCL256Architecture>(atm, config);