diff --git a/CMakeLists.txt b/CMakeLists.txt index 112e3b5..54d2362 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -58,6 +58,7 @@ install( MuonShield Magnet DecayVolume + SBT Trackers Calorimeter UpstreamTagger @@ -78,6 +79,7 @@ foreach( MuonShield Magnet DecayVolume + SBT Trackers Calorimeter UpstreamTagger diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 8cd3b93..6a39998 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -21,6 +21,7 @@ target_link_libraries( Target MuonShield DecayVolume + SBT UpstreamTagger Trackers Magnet diff --git a/src/SHiPGeometry.cpp b/src/SHiPGeometry.cpp index 7ff5d1a..f5b817e 100644 --- a/src/SHiPGeometry.cpp +++ b/src/SHiPGeometry.cpp @@ -15,6 +15,8 @@ #include "UpstreamTagger/SHiPUBTManager.h" #include "UpstreamTagger/UpstreamTaggerFactory.h" +#include "SBT/SBTFactory.h" + #include #include #include @@ -80,6 +82,18 @@ GeoPhysVol* SHiPGeometryBuilder::build() { world->add(new GeoTransform(decayVolumeTrf)); world->add(decayVolume); + // Build and place SBT (Surround Background Tagger). + // Wraps the decay volume — same Z range and centre (z = 32.92 to 83.32 m, + // centre 58.12 m). The intentional overlap with the DecayVolume is + // whitelisted in test_consistency. + SBTFactory sbtFactory(materials); + GeoPhysVol* sbt = sbtFactory.build(); + GeoTrf::Transform3D sbtTrf = GeoTrf::Translate3D(0.0, 0.0, 58.12 * 1000.0); + world->add(new GeoNameTag("/SHiP/sbt")); + world->add(new GeoIdentifierTag(9)); + world->add(new GeoTransform(sbtTrf)); + world->add(sbt); + // Build and place Trackers (container with 4 stations) // Container spans stations 1-4, centred appropriately TrackersFactory trackersFactory(materials); @@ -116,12 +130,13 @@ GeoPhysVol* SHiPGeometryBuilder::build() { world->add(new GeoTransform(timingDetectorTrf)); world->add(timingDetector); - // Build and place Calorimeter (container with ECAL front, ECAL back, HCAL) - // Full calorimeter spans Z: 96.87 to 99.77 m → centre: 98.32 m + // Build and place Calorimeter (ECAL + HCAL). + // The layer structure is driven by calo.cfg; the outer container dimensions + // and placement are fixed to match the SHiP subsystem envelope. CalorimeterFactory calorimeterFactory(materials); GeoPhysVol* calorimeter = calorimeterFactory.build(); GeoTrf::Transform3D calorimeterTrf = - GeoTrf::Translate3D(0.0, 0.0, 98.32 * 1000.0); // Convert m to mm + GeoTrf::Translate3D(0.0, 0.0, 98.32 * 1000.0); // 98.32 m, fixed envelope centre world->add(new GeoNameTag("/SHiP/calorimeter")); world->add(new GeoIdentifierTag(8)); world->add(new GeoTransform(calorimeterTrf)); diff --git a/src/SHiPMaterials.cpp b/src/SHiPMaterials.cpp index a644ea0..00c39d2 100644 --- a/src/SHiPMaterials.cpp +++ b/src/SHiPMaterials.cpp @@ -76,6 +76,10 @@ void SHiPMaterials::createElements() { // Aluminium m_elements["Aluminium"] = new GeoElement( "Aluminium", "Al", 13.0, 26.982 * GeoModelKernelUnits::g / GeoModelKernelUnits::mole); + + // Lead + m_elements["Lead"] = new GeoElement("Lead", "Pb", 82.0, + 207.2 * GeoModelKernelUnits::g / GeoModelKernelUnits::mole); } void SHiPMaterials::createMaterials() { @@ -179,6 +183,88 @@ void SHiPMaterials::createMaterials() { scintillator->add(m_elements["Hydrogen"], 0.085); scintillator->lock(); m_materials["Scintillator"] = scintillator; + + // Lead (density 11.34 g/cm³): pure Pb + GeoMaterial* lead = + new GeoMaterial("Lead", 11.34 * GeoModelKernelUnits::g / GeoModelKernelUnits::cm3); + lead->add(m_elements["Lead"], 1.0); + lead->lock(); + m_materials["Lead"] = lead; + + // PVT / polyvinyltoluene (density 1.032 g/cm³): C9H10 + // MW = 9*12.011 + 10*1.008 = 108.099 + 10.080 = 118.179 g/mol + { + const double awC = 12.011; + const double awH = 1.008; + const double mw = 9.0 * awC + 10.0 * awH; + GeoMaterial* pvt = + new GeoMaterial("PVT", 1.032 * GeoModelKernelUnits::g / GeoModelKernelUnits::cm3); + pvt->add(m_elements["Carbon"], 9.0 * awC / mw); + pvt->add(m_elements["Hydrogen"], 10.0 * awH / mw); + pvt->lock(); + m_materials["PVT"] = pvt; + } + + // Polystyrene (density 1.05 g/cm³): C8H8 + // MW = 8*12.011 + 8*1.008 = 96.088 + 8.064 = 104.152 g/mol + { + const double awC = 12.011; + const double awH = 1.008; + const double mw = 8.0 * awC + 8.0 * awH; + GeoMaterial* polystyrene = new GeoMaterial( + "Polystyrene", 1.05 * GeoModelKernelUnits::g / GeoModelKernelUnits::cm3); + polystyrene->add(m_elements["Carbon"], 8.0 * awC / mw); + polystyrene->add(m_elements["Hydrogen"], 8.0 * awH / mw); + polystyrene->lock(); + m_materials["Polystyrene"] = polystyrene; + } + + // Mylar / PET (density 1.39 g/cm³): C10H8O4, used for straw tube walls. + // MW = 10*12.011 + 8*1.008 + 4*15.999 = 120.11 + 8.064 + 63.996 = 192.170 g/mol + { + const double awC = 12.011; + const double awH = 1.008; + const double awO = 15.999; + const double mw = 10.0 * awC + 8.0 * awH + 4.0 * awO; + GeoMaterial* mylar = + new GeoMaterial("Mylar", 1.39 * GeoModelKernelUnits::g / GeoModelKernelUnits::cm3); + mylar->add(m_elements["Carbon"], 10.0 * awC / mw); + mylar->add(m_elements["Hydrogen"], 8.0 * awH / mw); + mylar->add(m_elements["Oxygen"], 4.0 * awO / mw); + mylar->lock(); + m_materials["Mylar"] = mylar; + } + + // Ar/CO2 70/30 by mass (density 1.56e-3 g/cm³), straw tube fill gas. + // CO2 mass-fraction is split into C and O by stoichiometry of the molecule. + { + const double awC = 12.011; + const double awO = 15.999; + const double mwCO2 = awC + 2.0 * awO; + const double fracAr = 0.70; + const double fracCO2 = 0.30; + GeoMaterial* arco2 = new GeoMaterial( + "ArCO2_70_30", 1.56e-3 * GeoModelKernelUnits::g / GeoModelKernelUnits::cm3); + arco2->add(m_elements["Argon"], fracAr); + arco2->add(m_elements["Carbon"], fracCO2 * awC / mwCO2); + arco2->add(m_elements["Oxygen"], fracCO2 * 2.0 * awO / mwCO2); + arco2->lock(); + m_materials["ArCO2_70_30"] = arco2; + } + + // LAB / Linear Alkylbenzene liquid scintillator (density 0.867 g/cm³). + // Representative formula C17.7H30.4 (n ≈ 12 in C6H5-CnH2n+1). + // Mass fractions: + // C: 17.7 * 12.011 / (17.7*12.011 + 30.4*1.008) = 212.59 / 243.23 = 0.8741 + // H: 30.4 * 1.008 / 243.23 = 0.1260 + { + GeoMaterial* lab = + new GeoMaterial("LAB", 0.867 * GeoModelKernelUnits::g / GeoModelKernelUnits::cm3); + lab->add(m_elements["Carbon"], 0.8741); + lab->add(m_elements["Hydrogen"], 0.1260); + lab->lock(); + m_materials["LAB"] = lab; + } } } // namespace SHiPGeometry diff --git a/subsystems/CMakeLists.txt b/subsystems/CMakeLists.txt index 00f0a7c..f1c47be 100644 --- a/subsystems/CMakeLists.txt +++ b/subsystems/CMakeLists.txt @@ -7,6 +7,7 @@ add_subdirectory(Target) add_subdirectory(MuonShield) add_subdirectory(Magnet) add_subdirectory(DecayVolume) +add_subdirectory(SBT) add_subdirectory(Trackers) add_subdirectory(Calorimeter) add_subdirectory(UpstreamTagger) diff --git a/subsystems/Calorimeter/CMakeLists.txt b/subsystems/Calorimeter/CMakeLists.txt index 89cceda..9eacd90 100644 --- a/subsystems/Calorimeter/CMakeLists.txt +++ b/subsystems/Calorimeter/CMakeLists.txt @@ -1,7 +1,26 @@ # SPDX-License-Identifier: LGPL-3.0-or-later # Copyright (C) CERN for the benefit of the SHiP Collaboration -add_library(Calorimeter src/CalorimeterFactory.cpp) +# Stage calo.cfg into the build directory so tests running from there find it. +configure_file( + ${CMAKE_CURRENT_SOURCE_DIR}/calo.cfg + ${CMAKE_CURRENT_BINARY_DIR}/calo.cfg + COPYONLY +) +# Also stage into the top-level build dir for test_builder and OverlapCheck. +configure_file( + ${CMAKE_CURRENT_SOURCE_DIR}/calo.cfg + ${CMAKE_BINARY_DIR}/calo.cfg + COPYONLY +) + +add_library( + Calorimeter + src/CalorimeterFactory.cpp + src/CalorimeterConfig.cpp + src/CaloBarLayer.cpp + src/CaloFibreHPLayer.cpp +) target_include_directories( Calorimeter @@ -13,6 +32,23 @@ target_include_directories( target_link_libraries(Calorimeter PUBLIC GeoModelCore::GeoModelKernel) +# Bake the absolute source-tree path as a compile-time fallback so the factory +# can always find calo.cfg even when CWD doesn't contain it. +target_compile_definitions( + Calorimeter + PRIVATE + # Source-tree fallback (always valid during development and CI builds). + CALO_CFG_DEFAULT_PATH="${CMAKE_CURRENT_SOURCE_DIR}/calo.cfg" + # Install-time fallback: set to the installed data directory so that + # deployed builds (where the source tree is absent) can find calo.cfg. + CALO_CFG_INSTALL_PATH="${CMAKE_INSTALL_PREFIX}/${CMAKE_INSTALL_DATADIR}/SHiPGeometry/calo.cfg" +) + +install( + FILES ${CMAKE_CURRENT_SOURCE_DIR}/calo.cfg + DESTINATION ${CMAKE_INSTALL_DATADIR}/SHiPGeometry +) + if(BUILD_TESTING) include(Catch) add_executable(test_calorimeter test_calorimeter.cpp) @@ -20,5 +56,7 @@ if(BUILD_TESTING) test_calorimeter PRIVATE Calorimeter SHiPGeometry Catch2::Catch2WithMain ) - catch_discover_tests(test_calorimeter) + catch_discover_tests(test_calorimeter + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + ) endif() diff --git a/subsystems/Calorimeter/calo.cfg b/subsystems/Calorimeter/calo.cfg new file mode 100644 index 0000000..648f631 --- /dev/null +++ b/subsystems/Calorimeter/calo.cfg @@ -0,0 +1,38 @@ +plate_xy_mm = 2160 +lead_thickness_mm = 3 +scint_thickness_mm = 10 +pvt_thickness_mm = 10 + +#How many calo modules in the x and y directions? +module_nx = 2 +module_ny = 3 + +module_pitch_x_mm = 0.0 +module_pitch_y_mm = 0.0 + +tol_x_mm = 10 +tol_y_mm = 10 +tol_z_mm = 10 + +# Layer codes: 1=WidePVT, 3=ThinPS, 5 HPL, 7=Lead, 8 split, maybe there is a better way to input these codes? +layers = 7,1,7,2,7,3,7,4,7,1,7,2,7,3,7,4,7,1,7,2,5,6,8,7,3,7,4,7,1,7,2,5,6,7,3,7,4,7,1,7,2,7,3,7,4,7,1,7,2,5,6,7,3,7,4,7,1,7,2,7,3,7,4,7,1,7,2,7,3,7,4,7,1,7,2,7,3,7,4,7,1,7,2,7,3,7,4 + +#HPLs +hpl_thickness_mm = 10 +fiber_diameter_mm = 1.2 +fiber_core_diameter_mm = 1.0 + + +airgap_mm = 1000 + + +#HCAL section + +gap_ecal_hcal_mm = 100; +#Layer codes are identical to above. 7 is iron though +layers2 = 7,1,7,2,7,1,7,2,7,1 +iron_thickness_mm = 170 + +detector_offset_z_mm = 96970.0 +detector_offset_x_mm = 0.0 +detector_offset_y_mm = 0.0 diff --git a/subsystems/Calorimeter/include/Calorimeter/CaloBarLayer.h b/subsystems/Calorimeter/include/Calorimeter/CaloBarLayer.h new file mode 100644 index 0000000..8dd3207 --- /dev/null +++ b/subsystems/Calorimeter/include/Calorimeter/CaloBarLayer.h @@ -0,0 +1,30 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) CERN for the benefit of the SHiP Collaboration + +#pragma once + +#include + +class GeoVPhysVol; +class GeoLogVol; + +namespace SHiPGeometry { + +/// Which transverse axis the bars are replicated along. +enum class BarAxis { AlongX, AlongY }; + +/** + * @brief Places an array of identical scintillator bars into a mother volume. + * + * Bars share a single GeoLogVol (reuse). They are spaced at @p pitch_mm + * centre-to-centre, centred on the mother origin in the transverse plane, + * and all placed at @p zCenter_mm along Z (local coordinates). + */ +class CaloBarLayer { + public: + static void place(GeoVPhysVol* mother, GeoLogVol* barLog, double pitch_mm, int nBars, + double zCenter_mm, const char* tagPrefix, int layerIndex, BarAxis axis, + const std::string& nameSuffix = ""); +}; + +} // namespace SHiPGeometry diff --git a/subsystems/Calorimeter/include/Calorimeter/CaloFibreHPLayer.h b/subsystems/Calorimeter/include/Calorimeter/CaloFibreHPLayer.h new file mode 100644 index 0000000..46a35b9 --- /dev/null +++ b/subsystems/Calorimeter/include/Calorimeter/CaloFibreHPLayer.h @@ -0,0 +1,28 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) CERN for the benefit of the SHiP Collaboration + +#pragma once + +#include + +class GeoVPhysVol; +class GeoMaterial; + +namespace SHiPGeometry { + +/** + * @brief Builds one Hadronic Pre-shower Layer (HPL) of scintillating fibres. + * + * The layer consists of an aluminium casing filled with three sublayers of + * cylindrical fibres (cladding + core). Fibres run along Y when + * @p fibresAlongY is true, along X otherwise. + */ +class CaloFibreHPLayer { + public: + static void build(GeoVPhysVol* mother, GeoMaterial* aluminiumMat, GeoMaterial* fibreMat, + const std::string& layerTag, double zCenter_mm, int layerIndex, + double casingXY_mm, double casingZ_mm, double fiberDiam_mm, + double fiberCoreDiam_mm, bool fibresAlongY, const std::string& nameSuffix); +}; + +} // namespace SHiPGeometry diff --git a/subsystems/Calorimeter/include/Calorimeter/CalorimeterConfig.h b/subsystems/Calorimeter/include/Calorimeter/CalorimeterConfig.h new file mode 100644 index 0000000..8803507 --- /dev/null +++ b/subsystems/Calorimeter/include/Calorimeter/CalorimeterConfig.h @@ -0,0 +1,77 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) CERN for the benefit of the SHiP Collaboration + +#pragma once + +#include +#include + +namespace SHiPGeometry { + +/** + * @brief Configuration for the SHiP calorimeter geometry. + * + * Populated by readCaloConfig() from a calo.cfg key=value file. + * + * Layer codes (layers / layers2 vectors): + * 1 = WidePVT bar layer, bars along X (H orientation) + * 2 = WidePVT bar layer, bars along Y (V orientation) + * 3 = ThinPS bar layer, bars along X (H orientation) + * 4 = ThinPS bar layer, bars along Y (V orientation) + * 5 = HPL fibre layer, fibres along Y + * 6 = HPL fibre layer, fibres along X + * 7 = absorber plate (Lead in ECAL section, Iron in HCAL section) + * 8 = air gap (no volume placed, just advances z cursor) + * + * layers → ECAL section + * layers2 → HCAL section (code 7 means Iron here) + */ +struct CalorimeterConfig { + // ECAL layer sequence + std::vector layers; + + // Individual layer thicknesses (mm) + double plate_xy_mm = 2160.0; + double lead_thickness_mm = 3.0; + double scint_thickness_mm = 10.0; + double pvt_thickness_mm = 10.0; + double hpl_thickness_mm = 10.0; + double fiber_diameter_mm = 1.2; + double fiber_core_diameter_mm = -1.0; // <0 → use fiber_diameter_mm + double airgap_mm = 1000.0; + + // HCAL layer sequence and absorber thickness + std::vector layers2; + double iron_thickness_mm = 170.0; + + // Gap between ECAL and HCAL stacks (mm) + double gap_ecal_hcal_mm = 0.0; + + // Module tiling + int module_nx = 1; + int module_ny = 1; + double module_pitch_x_mm = 0.0; // 0 → use plate_xy_mm + double module_pitch_y_mm = 0.0; // 0 → use plate_xy_mm + + // Tolerance (envelope padding, mm) + double tol_x_mm = 10.0; + double tol_y_mm = 10.0; + double tol_z_mm = 10.0; + + // Global position of the calorimeter stack centre in the world (mm) + double detector_offset_x_mm = 0.0; + double detector_offset_y_mm = 0.0; + double detector_offset_z_mm = 96970.0; + + // If true, the layer stack is centred at z=0 in the container volume. + bool center_stack = true; +}; + +/** + * @brief Parse a calo.cfg file and return a CalorimeterConfig. + * @throws std::runtime_error if the file cannot be opened or is missing + * mandatory fields. + */ +CalorimeterConfig readCaloConfig(const std::string& path); + +} // namespace SHiPGeometry diff --git a/subsystems/Calorimeter/include/Calorimeter/CalorimeterFactory.h b/subsystems/Calorimeter/include/Calorimeter/CalorimeterFactory.h index 6c42244..63d13bb 100644 --- a/subsystems/Calorimeter/include/Calorimeter/CalorimeterFactory.h +++ b/subsystems/Calorimeter/include/Calorimeter/CalorimeterFactory.h @@ -3,62 +3,61 @@ #pragma once +#include + class GeoPhysVol; namespace SHiPGeometry { class SHiPMaterials; +struct CalorimeterConfig; /** - * @brief Factory for the Calorimeter (electromagnetic + hadronic calorimeter) geometry + * @brief Factory for the Calorimeter (ECAL + HCAL) geometry. + * + * Creates a fixed-size container volume matching the SHiP envelope + * (3.00 × 3.50 × 1.45 m half-sizes, centred at Z = 98 320 mm) and fills + * it with the real layer-by-layer geometry driven by calo.cfg. * - * Creates a container with three sections: - * - ECAL front: Z 96.87-97.07 m → centre 96.97 m, half-length 0.10 m, 2.25×3.50 m - * - ECAL back: Z 98.07-98.67 m → centre 98.37 m, half-length 0.30 m, 2.75×3.50 m - * - HCAL: Z 98.77-99.77 m → centre 99.27 m, half-length 0.50 m, 3.00×3.50 m + * The config file is resolved at build() time: + * 1. "calo.cfg" relative to the current working directory (works when + * running from the build directory, where CMake stages the file). + * 2. The absolute source-tree path baked in at compile time via + * CALO_CFG_DEFAULT_PATH (always valid during development). */ class CalorimeterFactory { public: - explicit CalorimeterFactory(SHiPMaterials& materials); + explicit CalorimeterFactory(SHiPMaterials& materials, std::string configPath = "calo.cfg"); ~CalorimeterFactory() = default; + /** Build and return the calorimeter container volume. */ + GeoPhysVol* build(); + /** - * @brief Build the Calorimeter geometry - * @return Pointer to container volume with ECAL front, ECAL back, and HCAL + * @brief Compute the total Z extent of one ECAL+gap+HCAL stack (mm). + * + * Public so tests and placement code can query the stack thickness + * independently of a full build(). */ - GeoPhysVol* build(); + static double totalStackZ(const CalorimeterConfig& cfg); + + /** Return the config path that will actually be opened (after resolution). */ + std::string resolvedConfigPath() const; private: SHiPMaterials& m_materials; - - // Component dimensions (convert m to mm) - // Note: GeoModel uses mm internally, so 1m = 1000mm - - // ECAL front - static constexpr double s_ecalFrontHalfX = 2250.0; // 2.25 m - static constexpr double s_ecalFrontHalfY = 3500.0; // 3.50 m - static constexpr double s_ecalFrontHalfZ = 100.0; // 0.10 m - static constexpr double s_ecalFrontZ = 96970.0; // 96.97 m - - // ECAL back - static constexpr double s_ecalBackHalfX = 2750.0; // 2.75 m - static constexpr double s_ecalBackHalfY = 3500.0; // 3.50 m - static constexpr double s_ecalBackHalfZ = 300.0; // 0.30 m - static constexpr double s_ecalBackZ = 98370.0; // 98.37 m - - // HCAL - static constexpr double s_hcalHalfX = 3000.0; // 3.00 m - static constexpr double s_hcalHalfY = 3500.0; // 3.50 m - static constexpr double s_hcalHalfZ = 500.0; // 0.50 m - static constexpr double s_hcalZ = 99270.0; // 99.27 m - - // Container dimensions (spans all components) - // From CSV: ECAL front 96.87-97.07, ECAL back 98.07-98.67, HCAL 98.77-99.77 - // Full span: 96.87 to 99.77 m → centre: 98.32 m, half-length: 1.45 m - static constexpr double s_containerHalfX = s_hcalHalfX; // Use largest X - static constexpr double s_containerHalfY = s_hcalHalfY; // Use largest Y - static constexpr double s_containerHalfZ = 1450.0; // 1.45 m - static constexpr double s_containerCentreZ = 98320.0; // 98.32 m + std::string m_configPath; + + /** Place one NX×NY tiled stack of layers inside @p container. */ + void buildStack(GeoPhysVol* container, const CalorimeterConfig& cfg, int mx, int my, + double offsetX, double offsetY) const; + + // ── Fixed container dimensions (mm) ───────────────────────────────── + // These match the SHiP subsystem envelope from subsystem_envelopes.csv + // and must not change — tests and the consistency check depend on them. + static constexpr double s_containerHalfX = 3000.0; // 3.00 m + static constexpr double s_containerHalfY = 3500.0; // 3.50 m + static constexpr double s_containerHalfZ = 1450.0; // 1.45 m }; } // namespace SHiPGeometry diff --git a/subsystems/Calorimeter/src/CaloBarLayer.cpp b/subsystems/Calorimeter/src/CaloBarLayer.cpp new file mode 100644 index 0000000..fb49757 --- /dev/null +++ b/subsystems/Calorimeter/src/CaloBarLayer.cpp @@ -0,0 +1,43 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) CERN for the benefit of the SHiP Collaboration + +#include "Calorimeter/CaloBarLayer.h" + +#include +#include +#include +#include +#include + +#include + +namespace SHiPGeometry { + +using namespace GeoModelKernelUnits; + +void CaloBarLayer::place(GeoVPhysVol* mother, GeoLogVol* barLog, double pitch_mm, int nBars, + double zCenter_mm, const char* tagPrefix, int layerIndex, BarAxis axis, + const std::string& nameSuffix) { + const double pitch = pitch_mm * mm; + const double s0 = -0.5 * (nBars - 1) * pitch; + + for (int i = 0; i < nBars; ++i) { + const double s = s0 + i * pitch; + + double x = 0.0, y = 0.0; + if (axis == BarAxis::AlongX) + x = s; + else + y = s; + + const std::string name = std::string(tagPrefix) + "_L" + std::to_string(layerIndex) + "_B" + + std::to_string(i) + nameSuffix; + + auto* barPhys = new GeoPhysVol(barLog); + mother->add(new GeoNameTag(name.c_str())); + mother->add(new GeoTransform(GeoTrf::Translate3D(x, y, zCenter_mm * mm))); + mother->add(barPhys); + } +} + +} // namespace SHiPGeometry diff --git a/subsystems/Calorimeter/src/CaloFibreHPLayer.cpp b/subsystems/Calorimeter/src/CaloFibreHPLayer.cpp new file mode 100644 index 0000000..a89024b --- /dev/null +++ b/subsystems/Calorimeter/src/CaloFibreHPLayer.cpp @@ -0,0 +1,97 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) CERN for the benefit of the SHiP Collaboration + +#include "Calorimeter/CaloFibreHPLayer.h" + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +namespace SHiPGeometry { + +using namespace GeoModelKernelUnits; + +void CaloFibreHPLayer::build(GeoVPhysVol* mother, GeoMaterial* aluminiumMat, GeoMaterial* fibreMat, + const std::string& layerTag, double zCenter_mm, int layerIndex, + double casingXY_mm, double casingZ_mm, double fiberDiam_mm, + double fiberCoreDiam_mm, bool fibresAlongY, + const std::string& nameSuffix) { + const double casingXY = casingXY_mm * mm; + const double casingZ = casingZ_mm * mm; + + // ── aluminium casing ────────────────────────────────────────────────── + auto* casingShape = new GeoBox(0.5 * casingXY, 0.5 * casingXY, 0.5 * casingZ); + auto* casingLog = new GeoLogVol("HPL_CasingLog", casingShape, aluminiumMat); + auto* casingPhys = new GeoPhysVol(casingLog); + + mother->add(new GeoNameTag((layerTag + "_HPL_Casing" + nameSuffix).c_str())); + mother->add(new GeoTransform(GeoTrf::Translate3D(0.0, 0.0, zCenter_mm * mm))); + mother->add(casingPhys); + + // ── fibre geometry ──────────────────────────────────────────────────── + const double rOuter = 0.5 * fiberDiam_mm * mm; + const double rCore = 0.5 * fiberCoreDiam_mm * mm; + + if (rCore <= 0.0 || rCore > rOuter) + throw std::runtime_error( + "CaloFibreHPLayer: invalid core diameter " + "(must be >0 and <= outer diameter)"); + + // GeoTube axis is Z; rotate so fibres run along Y or X + const double halfLen = 0.5 * casingXY; + GeoTrf::Transform3D rotAxis = + fibresAlongY ? GeoTrf::Transform3D(GeoTrf::RotateX3D(90.0 * deg)) // Z → Y + : GeoTrf::Transform3D(GeoTrf::RotateY3D(-90.0 * deg)); // Z → X + + auto* cladShape = new GeoTube(0.0, rOuter, halfLen); + auto* coreShape = new GeoTube(0.0, rCore, halfLen); + auto* cladLog = new GeoLogVol("HPL_CladLog", cladShape, fibreMat); + auto* coreLog = new GeoLogVol("HPL_FiberCoreLog", coreShape, fibreMat); + + // ── fibre tiling: three staggered sublayers ─────────────────────────── + const double pitch = 2.0 * rOuter; // tight packing + const double dxMax = 0.5 * pitch; // middle-sublayer shift + const double usable = casingXY - 2.0 * rOuter - dxMax; + const int nFib = std::max(1, static_cast(std::floor(usable / pitch)) + 1); + const double x0 = -0.5 * ((nFib - 1) * pitch + dxMax); + + const double dx[3] = {0.0, 0.5 * pitch, 0.0}; + const double dz = std::max(0.0, 0.5 * casingZ - rOuter); + const double zLocal[3] = {-dz, 0.0, +dz}; + + const std::string orient = fibresAlongY ? "V_L" : "H_L"; + + for (int s = 0; s < 3; ++s) { + for (int i = 0; i < nFib; ++i) { + const double pack = x0 + i * pitch + dx[s]; + const double zl = zLocal[s]; + + const std::string baseName = layerTag + "_HPL_" + orient + std::to_string(layerIndex) + + "_S" + std::to_string(s) + "_F" + std::to_string(i) + + nameSuffix; + + // cladding physical volume; core is a child of it + auto* cladPhys = new GeoPhysVol(cladLog); + cladPhys->add(new GeoNameTag((baseName + "_Core").c_str())); + cladPhys->add(new GeoPhysVol(coreLog)); + + const double xPos = fibresAlongY ? pack : 0.0; + const double yPos = fibresAlongY ? 0.0 : pack; + + casingPhys->add(new GeoNameTag((baseName + "_Clad").c_str())); + casingPhys->add(new GeoTransform(GeoTrf::Translate3D(xPos, yPos, zl) * rotAxis)); + casingPhys->add(cladPhys); + } + } +} + +} // namespace SHiPGeometry diff --git a/subsystems/Calorimeter/src/CalorimeterConfig.cpp b/subsystems/Calorimeter/src/CalorimeterConfig.cpp new file mode 100644 index 0000000..fdf04dd --- /dev/null +++ b/subsystems/Calorimeter/src/CalorimeterConfig.cpp @@ -0,0 +1,144 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) CERN for the benefit of the SHiP Collaboration + +#include "Calorimeter/CalorimeterConfig.h" + +#include +#include +#include +#include +#include +#include + +namespace SHiPGeometry { + +namespace { + +std::string trim(std::string s) { + auto notSpace = [](unsigned char c) { return !std::isspace(c); }; + s.erase(s.begin(), std::find_if(s.begin(), s.end(), notSpace)); + s.erase(std::find_if(s.rbegin(), s.rend(), notSpace).base(), s.end()); + return s; +} + +std::vector parseIntList(const std::string& s) { + std::vector out; + std::stringstream ss(s); + std::string token; + while (std::getline(ss, token, ',')) { + token = trim(token); + // strip trailing semicolons (calo.cfg has e.g. "gap_ecal_hcal_mm = 100;") + if (!token.empty() && token.back() == ';') + token.pop_back(); + if (!token.empty()) + out.push_back(std::stoi(token)); + } + return out; +} + +double toDouble(std::string v) { + // strip optional trailing semicolon + if (!v.empty() && v.back() == ';') + v.pop_back(); + return std::stod(v); +} + +int toInt(std::string v) { + if (!v.empty() && v.back() == ';') + v.pop_back(); + return std::stoi(v); +} + +} // namespace + +CalorimeterConfig readCaloConfig(const std::string& path) { + CalorimeterConfig cfg; + + std::ifstream in(path); + if (!in) + throw std::runtime_error("CalorimeterConfig: cannot open: " + path); + + std::string line; + while (std::getline(in, line)) { + // strip comments + auto hash = line.find('#'); + if (hash != std::string::npos) + line = line.substr(0, hash); + line = trim(line); + if (line.empty()) + continue; + + auto eq = line.find('='); + if (eq == std::string::npos) + continue; + + auto key = trim(line.substr(0, eq)); + auto val = trim(line.substr(eq + 1)); + + if (key == "layers") + cfg.layers = parseIntList(val); + else if (key == "layers2") + cfg.layers2 = parseIntList(val); + else if (key == "plate_xy_mm") + cfg.plate_xy_mm = toDouble(val); + else if (key == "lead_thickness_mm") + cfg.lead_thickness_mm = toDouble(val); + else if (key == "scint_thickness_mm") + cfg.scint_thickness_mm = toDouble(val); + else if (key == "pvt_thickness_mm") + cfg.pvt_thickness_mm = toDouble(val); + else if (key == "hpl_thickness_mm") + cfg.hpl_thickness_mm = toDouble(val); + else if (key == "fiber_diameter_mm") + cfg.fiber_diameter_mm = toDouble(val); + else if (key == "fiber_core_diameter_mm") + cfg.fiber_core_diameter_mm = toDouble(val); + else if (key == "airgap_mm") + cfg.airgap_mm = toDouble(val); + else if (key == "iron_thickness_mm") + cfg.iron_thickness_mm = toDouble(val); + else if (key == "gap_ecal_hcal_mm") + cfg.gap_ecal_hcal_mm = toDouble(val); + else if (key == "module_nx") + cfg.module_nx = toInt(val); + else if (key == "module_ny") + cfg.module_ny = toInt(val); + else if (key == "module_pitch_x_mm") + cfg.module_pitch_x_mm = toDouble(val); + else if (key == "module_pitch_y_mm") + cfg.module_pitch_y_mm = toDouble(val); + else if (key == "tol_x_mm") + cfg.tol_x_mm = toDouble(val); + else if (key == "tol_y_mm") + cfg.tol_y_mm = toDouble(val); + else if (key == "tol_z_mm") + cfg.tol_z_mm = toDouble(val); + else if (key == "detector_offset_x_mm") + cfg.detector_offset_x_mm = toDouble(val); + else if (key == "detector_offset_y_mm") + cfg.detector_offset_y_mm = toDouble(val); + else if (key == "detector_offset_z_mm") + cfg.detector_offset_z_mm = toDouble(val); + else if (key == "center_stack") { + std::string v = val; + std::transform(v.begin(), v.end(), v.begin(), ::tolower); + cfg.center_stack = (v == "1" || v == "true" || v == "yes" || v == "on"); + } + // unknown keys are silently ignored to allow forward-compatibility + } + + if (cfg.layers.empty()) + throw std::runtime_error("CalorimeterConfig: 'layers' must be defined in " + path); + + // Default fibre-core diameter to outer diameter when not set + if (cfg.fiber_core_diameter_mm <= 0) + cfg.fiber_core_diameter_mm = cfg.fiber_diameter_mm; + + if (cfg.fiber_core_diameter_mm > cfg.fiber_diameter_mm) + throw std::runtime_error( + "CalorimeterConfig: fiber_core_diameter_mm cannot exceed fiber_diameter_mm"); + + return cfg; +} + +} // namespace SHiPGeometry diff --git a/subsystems/Calorimeter/src/CalorimeterFactory.cpp b/subsystems/Calorimeter/src/CalorimeterFactory.cpp index ef3c9d3..8e53257 100644 --- a/subsystems/Calorimeter/src/CalorimeterFactory.cpp +++ b/subsystems/Calorimeter/src/CalorimeterFactory.cpp @@ -3,6 +3,9 @@ #include "Calorimeter/CalorimeterFactory.h" +#include "Calorimeter/CaloBarLayer.h" +#include "Calorimeter/CaloFibreHPLayer.h" +#include "Calorimeter/CalorimeterConfig.h" #include "SHiPGeometry/SHiPMaterials.h" #include @@ -12,59 +15,382 @@ #include #include #include +#include + +#include +#include +#include + +// Absolute fallback path baked in by CMake so out-of-source builds always +// find calo.cfg even when the CWD doesn't contain a copy of it. +#ifndef CALO_CFG_DEFAULT_PATH +#define CALO_CFG_DEFAULT_PATH "calo.cfg" +#endif +// Install-time data directory path, set by CMake during install configuration. +#ifndef CALO_CFG_INSTALL_PATH +#define CALO_CFG_INSTALL_PATH "" +#endif namespace SHiPGeometry { -CalorimeterFactory::CalorimeterFactory(SHiPMaterials& materials) : m_materials(materials) {} +using namespace GeoModelKernelUnits; + +// ── file-scope helper ──────────────────────────────────────────────────────── + +static std::string resolveCfgPath(const std::string& path) { + if (!path.empty() && path[0] == '/') + return path; // already absolute + if (std::ifstream(path).good()) + return path; // found relative to CWD + const std::string srcFallback = CALO_CFG_DEFAULT_PATH; + if (std::ifstream(srcFallback).good()) + return srcFallback; // source-tree fallback + const std::string installFallback = CALO_CFG_INSTALL_PATH; + if (!installFallback.empty() && std::ifstream(installFallback).good()) + return installFallback; // installed data dir + return path; // give up — readCaloConfig will emit the error +} + +// ── constructor / accessors ────────────────────────────────────────────────── + +CalorimeterFactory::CalorimeterFactory(SHiPMaterials& materials, std::string configPath) + : m_materials(materials), m_configPath(std::move(configPath)) {} + +std::string CalorimeterFactory::resolvedConfigPath() const { + return resolveCfgPath(m_configPath); +} + +// ── totalStackZ ────────────────────────────────────────────────────────────── + +double CalorimeterFactory::totalStackZ(const CalorimeterConfig& cfg) { + double total = 0.0; + for (int code : cfg.layers) { + if (code == 7) + total += cfg.lead_thickness_mm; + else if (code >= 1 && code <= 4) + total += cfg.scint_thickness_mm; + else if (code == 5 || code == 6) + total += cfg.hpl_thickness_mm; + else if (code == 8) + total += cfg.airgap_mm; + else + throw std::runtime_error("CalorimeterFactory: unknown layer code in 'layers': " + + std::to_string(code)); + } + total += cfg.gap_ecal_hcal_mm; + for (int code : cfg.layers2) { + if (code == 7) + total += cfg.iron_thickness_mm; + else if (code >= 1 && code <= 4) + total += cfg.scint_thickness_mm; + else if (code == 5 || code == 6) + total += cfg.hpl_thickness_mm; + else if (code == 8) + total += cfg.airgap_mm; + else + throw std::runtime_error("CalorimeterFactory: unknown layer code in 'layers2': " + + std::to_string(code)); + } + return total; +} + +// ── build ──────────────────────────────────────────────────────────────────── GeoPhysVol* CalorimeterFactory::build() { - const GeoMaterial* air = m_materials.requireMaterial("Air"); + const CalorimeterConfig cfg = readCaloConfig(resolveCfgPath(m_configPath)); + + GeoMaterial* air = m_materials.requireMaterial("Air"); - // Create container volume that spans all calorimeter components + // Fixed-size container — must match the SHiP subsystem envelope so that + // the geometry consistency tests and the overlap check pass. auto* containerBox = new GeoBox(s_containerHalfX, s_containerHalfY, s_containerHalfZ); auto* containerLog = new GeoLogVol("/SHiP/calorimeter", containerBox, air); auto* containerPhys = new GeoPhysVol(containerLog); - // Create and place ECAL front - auto* ecalFrontBox = new GeoBox(s_ecalFrontHalfX, s_ecalFrontHalfY, s_ecalFrontHalfZ); - auto* ecalFrontLog = new GeoLogVol("/SHiP/calorimeter/ecal_front", ecalFrontBox, air); - auto* ecalFrontPhys = new GeoPhysVol(ecalFrontLog); + // Module pitch (0 → use plate size, modules touch) + const double pitchX = cfg.module_pitch_x_mm > 0.0 ? cfg.module_pitch_x_mm : cfg.plate_xy_mm; + const double pitchY = cfg.module_pitch_y_mm > 0.0 ? cfg.module_pitch_y_mm : cfg.plate_xy_mm; + + // ── Validate that the config fits inside the fixed container ───────────── + const double stackZ = totalStackZ(cfg); + const double maxHalfX = 0.5 * cfg.plate_xy_mm + 0.5 * (cfg.module_nx - 1) * pitchX; + const double maxHalfY = 0.5 * cfg.plate_xy_mm + 0.5 * (cfg.module_ny - 1) * pitchY; + if (stackZ > 2.0 * s_containerHalfZ) + throw std::runtime_error("CalorimeterFactory: calo.cfg total stack Z (" + + std::to_string(stackZ) + " mm) exceeds container half-Z*2 (" + + std::to_string(2.0 * s_containerHalfZ) + + " mm). Reduce layer thicknesses or number of layers."); + if (maxHalfX > s_containerHalfX) + throw std::runtime_error("CalorimeterFactory: calo.cfg module array half-X (" + + std::to_string(maxHalfX) + " mm) exceeds container half-X (" + + std::to_string(s_containerHalfX) + + " mm). Reduce plate_xy_mm, module_nx, or module_pitch_x_mm."); + if (maxHalfY > s_containerHalfY) + throw std::runtime_error("CalorimeterFactory: calo.cfg module array half-Y (" + + std::to_string(maxHalfY) + " mm) exceeds container half-Y (" + + std::to_string(s_containerHalfY) + + " mm). Reduce plate_xy_mm, module_ny, or module_pitch_y_mm."); + + const double x0 = -0.5 * (cfg.module_nx - 1) * pitchX; + const double y0 = -0.5 * (cfg.module_ny - 1) * pitchY; - double ecalFrontRelativeZ = s_ecalFrontZ - s_containerCentreZ; - GeoTrf::Transform3D ecalFrontTrf = GeoTrf::Translate3D(0.0, 0.0, ecalFrontRelativeZ); + for (int iy = 0; iy < cfg.module_ny; ++iy) + for (int ix = 0; ix < cfg.module_nx; ++ix) + buildStack(containerPhys, cfg, ix, iy, x0 + ix * pitchX, y0 + iy * pitchY); - containerPhys->add(new GeoNameTag("/SHiP/calorimeter/ecal_front")); - containerPhys->add(new GeoIdentifierTag(0)); - containerPhys->add(new GeoTransform(ecalFrontTrf)); - containerPhys->add(ecalFrontPhys); + return containerPhys; +} - // Create and place ECAL back - auto* ecalBackBox = new GeoBox(s_ecalBackHalfX, s_ecalBackHalfY, s_ecalBackHalfZ); - auto* ecalBackLog = new GeoLogVol("/SHiP/calorimeter/ecal_back", ecalBackBox, air); - auto* ecalBackPhys = new GeoPhysVol(ecalBackLog); +// ── buildStack ─────────────────────────────────────────────────────────────── - double ecalBackRelativeZ = s_ecalBackZ - s_containerCentreZ; - GeoTrf::Transform3D ecalBackTrf = GeoTrf::Translate3D(0.0, 0.0, ecalBackRelativeZ); +void CalorimeterFactory::buildStack(GeoPhysVol* container, const CalorimeterConfig& cfg, int mx, + int my, double offX, double offY) const { + GeoMaterial* leadMat = m_materials.requireMaterial("Lead"); + GeoMaterial* ironMat = m_materials.requireMaterial("Iron"); + GeoMaterial* pvtMat = m_materials.requireMaterial("PVT"); + GeoMaterial* psMat = m_materials.requireMaterial("Polystyrene"); + GeoMaterial* alMat = m_materials.requireMaterial("Aluminium"); + GeoMaterial* airMat = m_materials.requireMaterial("Air"); - containerPhys->add(new GeoNameTag("/SHiP/calorimeter/ecal_back")); - containerPhys->add(new GeoIdentifierTag(1)); - containerPhys->add(new GeoTransform(ecalBackTrf)); - containerPhys->add(ecalBackPhys); + const double plateXY = cfg.plate_xy_mm * mm; + const double leadZ = cfg.lead_thickness_mm * mm; + const double scintZ = cfg.scint_thickness_mm * mm; + const double hplZ = cfg.hpl_thickness_mm * mm; + const double ironZ = cfg.iron_thickness_mm * mm; + const double airGapZ = cfg.airgap_mm * mm; - // Create and place HCAL - auto* hcalBox = new GeoBox(s_hcalHalfX, s_hcalHalfY, s_hcalHalfZ); - auto* hcalLog = new GeoLogVol("/SHiP/calorimeter/hcal", hcalBox, air); - auto* hcalPhys = new GeoPhysVol(hcalLog); + const double wideW = 60.0 * mm; + const double thinW = 10.0 * mm; - double hcalRelativeZ = s_hcalZ - s_containerCentreZ; - GeoTrf::Transform3D hcalTrf = GeoTrf::Translate3D(0.0, 0.0, hcalRelativeZ); + const std::string mtag = "_MX" + std::to_string(mx) + "Y" + std::to_string(my); - containerPhys->add(new GeoNameTag("/SHiP/calorimeter/hcal")); - containerPhys->add(new GeoIdentifierTag(2)); - containerPhys->add(new GeoTransform(hcalTrf)); - containerPhys->add(hcalPhys); + // Reusable LogVols — GeoModel shares them across GeoPhysVol instances + auto* leadLog = new GeoLogVol("/SHiP/calorimeter/lead_plate" + mtag, + new GeoBox(0.5 * plateXY, 0.5 * plateXY, 0.5 * leadZ), leadMat); + auto* ironLog = new GeoLogVol("/SHiP/calorimeter/iron_plate" + mtag, + new GeoBox(0.5 * plateXY, 0.5 * plateXY, 0.5 * ironZ), ironMat); + auto* wideHLog = new GeoLogVol("/SHiP/calorimeter/wide_pvt_h" + mtag, + new GeoBox(0.5 * plateXY, 0.5 * wideW, 0.5 * scintZ), pvtMat); + auto* wideVLog = new GeoLogVol("/SHiP/calorimeter/wide_pvt_v" + mtag, + new GeoBox(0.5 * wideW, 0.5 * plateXY, 0.5 * scintZ), pvtMat); + auto* thinHLog = new GeoLogVol("/SHiP/calorimeter/thin_ps_h" + mtag, + new GeoBox(0.5 * plateXY, 0.5 * thinW, 0.5 * scintZ), psMat); + auto* thinVLog = new GeoLogVol("/SHiP/calorimeter/thin_ps_v" + mtag, + new GeoBox(0.5 * thinW, 0.5 * plateXY, 0.5 * scintZ), psMat); - return containerPhys; + // z cursor: start at -halfContainerZ if centering, else 0 + // We centre within the fixed container (halfZ = s_containerHalfZ). + const double totalZ = totalStackZ(cfg); + double zCursor = cfg.center_stack ? -0.5 * totalZ * mm : -s_containerHalfZ; + + int layerId = 0; // monotonic identifier used for GeoIdentifierTag + + // Lambda: wrap a layer in a named air-box envelope inside the container + auto makeEnv = [&](const std::string& name, double halfZ, double zCenter) -> GeoPhysVol* { + auto* s = new GeoBox(0.5 * plateXY, 0.5 * plateXY, halfZ); + auto* l = new GeoLogVol(name + "_env", s, airMat); + auto* p = new GeoPhysVol(l); + container->add(new GeoNameTag(name.c_str())); + container->add(new GeoIdentifierTag(layerId++)); + container->add(new GeoTransform(GeoTrf::Translate3D(offX, offY, zCenter))); + container->add(p); + return p; + }; + + int iWideH = 0, iWideV = 0, iThinH = 0, iThinV = 0, iHPL = 0; + int gl = 0, sl = 0; + + // ══════════════════════════════════════════ + // ECAL section (code 7 → Lead) + // ══════════════════════════════════════════ + for (int code : cfg.layers) { + if (code == 7) { + const std::string n = "/SHiP/calorimeter/ecal/gl" + std::to_string(gl) + "_lead" + mtag; + auto* env = makeEnv(n, 0.5 * leadZ, zCursor + 0.5 * leadZ); + env->add(new GeoNameTag(n.c_str())); + env->add(new GeoIdentifierTag(0)); + env->add(new GeoTransform(GeoTrf::Translate3D(0, 0, 0))); + env->add(new GeoPhysVol(leadLog)); + zCursor += leadZ; + ++gl; + + } else if (code == 1) { + const std::string n = "/SHiP/calorimeter/ecal/gl" + std::to_string(gl) + "_sl" + + std::to_string(sl) + "_wide_pvt_h" + mtag; + auto* env = makeEnv(n, 0.5 * scintZ, zCursor + 0.5 * scintZ); + CaloBarLayer::place(env, wideHLog, 60.0, 36, 0.0, n.c_str(), iWideH, BarAxis::AlongY, + mtag); + zCursor += scintZ; + ++iWideH; + ++gl; + ++sl; + + } else if (code == 2) { + const std::string n = "/SHiP/calorimeter/ecal/gl" + std::to_string(gl) + "_sl" + + std::to_string(sl) + "_wide_pvt_v" + mtag; + auto* env = makeEnv(n, 0.5 * scintZ, zCursor + 0.5 * scintZ); + CaloBarLayer::place(env, wideVLog, 60.0, 36, 0.0, n.c_str(), iWideV, BarAxis::AlongX, + mtag); + zCursor += scintZ; + ++iWideV; + ++gl; + ++sl; + + } else if (code == 3) { + const std::string n = "/SHiP/calorimeter/ecal/gl" + std::to_string(gl) + "_sl" + + std::to_string(sl) + "_thin_ps_h" + mtag; + auto* env = makeEnv(n, 0.5 * scintZ, zCursor + 0.5 * scintZ); + CaloBarLayer::place(env, thinHLog, 10.0, 216, 0.0, n.c_str(), iThinH, BarAxis::AlongY, + mtag); + zCursor += scintZ; + ++iThinH; + ++gl; + ++sl; + + } else if (code == 4) { + const std::string n = "/SHiP/calorimeter/ecal/gl" + std::to_string(gl) + "_sl" + + std::to_string(sl) + "_thin_ps_v" + mtag; + auto* env = makeEnv(n, 0.5 * scintZ, zCursor + 0.5 * scintZ); + CaloBarLayer::place(env, thinVLog, 10.0, 216, 0.0, n.c_str(), iThinV, BarAxis::AlongX, + mtag); + zCursor += scintZ; + ++iThinV; + ++gl; + ++sl; + + } else if (code == 5) { + const std::string n = "/SHiP/calorimeter/ecal/gl" + std::to_string(gl) + "_sl" + + std::to_string(sl) + "_hpl_y" + mtag; + auto* env = makeEnv(n, 0.5 * hplZ, zCursor + 0.5 * hplZ); + CaloFibreHPLayer::build(env, alMat, psMat, n, 0.0, iHPL, cfg.plate_xy_mm, + cfg.hpl_thickness_mm, cfg.fiber_diameter_mm, + cfg.fiber_core_diameter_mm, true, mtag); + zCursor += hplZ; + ++iHPL; + ++gl; + ++sl; + + } else if (code == 6) { + const std::string n = "/SHiP/calorimeter/ecal/gl" + std::to_string(gl) + "_sl" + + std::to_string(sl) + "_hpl_x" + mtag; + auto* env = makeEnv(n, 0.5 * hplZ, zCursor + 0.5 * hplZ); + CaloFibreHPLayer::build(env, alMat, psMat, n, 0.0, iHPL, cfg.plate_xy_mm, + cfg.hpl_thickness_mm, cfg.fiber_diameter_mm, + cfg.fiber_core_diameter_mm, false, mtag); + zCursor += hplZ; + ++iHPL; + ++gl; + ++sl; + + } else if (code == 8) { + zCursor += airGapZ; + + } else { + throw std::runtime_error("CalorimeterFactory: unknown layer code in 'layers': " + + std::to_string(code)); + } + } + + // ECAL–HCAL gap + zCursor += cfg.gap_ecal_hcal_mm * mm; + + // ══════════════════════════════════════════ + // HCAL section (code 7 → Iron) + // ══════════════════════════════════════════ + int iIron = 0; + gl = 0; + sl = 0; + + for (int code : cfg.layers2) { + if (code == 7) { + const std::string tag = "/SHiP/calorimeter/hcal/gl" + std::to_string(gl) + "_iron_" + + std::to_string(iIron) + mtag; + container->add(new GeoNameTag(tag.c_str())); + container->add(new GeoIdentifierTag(layerId++)); + container->add( + new GeoTransform(GeoTrf::Translate3D(offX, offY, zCursor + 0.5 * ironZ))); + container->add(new GeoPhysVol(ironLog)); + zCursor += ironZ; + ++iIron; + ++gl; + + } else if (code == 1) { + const std::string n = "/SHiP/calorimeter/hcal/gl" + std::to_string(gl) + "_sl" + + std::to_string(sl) + "_wide_pvt_h" + mtag; + auto* env = makeEnv(n, 0.5 * scintZ, zCursor + 0.5 * scintZ); + CaloBarLayer::place(env, wideHLog, 60.0, 36, 0.0, n.c_str(), iWideH, BarAxis::AlongY, + mtag); + zCursor += scintZ; + ++iWideH; + ++gl; + ++sl; + + } else if (code == 2) { + const std::string n = "/SHiP/calorimeter/hcal/gl" + std::to_string(gl) + "_sl" + + std::to_string(sl) + "_wide_pvt_v" + mtag; + auto* env = makeEnv(n, 0.5 * scintZ, zCursor + 0.5 * scintZ); + CaloBarLayer::place(env, wideVLog, 60.0, 36, 0.0, n.c_str(), iWideV, BarAxis::AlongX, + mtag); + zCursor += scintZ; + ++iWideV; + ++gl; + ++sl; + + } else if (code == 3) { + const std::string n = "/SHiP/calorimeter/hcal/gl" + std::to_string(gl) + "_sl" + + std::to_string(sl) + "_thin_ps_h" + mtag; + auto* env = makeEnv(n, 0.5 * scintZ, zCursor + 0.5 * scintZ); + CaloBarLayer::place(env, thinHLog, 10.0, 216, 0.0, n.c_str(), iThinH, BarAxis::AlongY, + mtag); + zCursor += scintZ; + ++iThinH; + ++gl; + ++sl; + + } else if (code == 4) { + const std::string n = "/SHiP/calorimeter/hcal/gl" + std::to_string(gl) + "_sl" + + std::to_string(sl) + "_thin_ps_v" + mtag; + auto* env = makeEnv(n, 0.5 * scintZ, zCursor + 0.5 * scintZ); + CaloBarLayer::place(env, thinVLog, 10.0, 216, 0.0, n.c_str(), iThinV, BarAxis::AlongX, + mtag); + zCursor += scintZ; + ++iThinV; + ++gl; + ++sl; + + } else if (code == 5) { + const std::string n = "/SHiP/calorimeter/hcal/gl" + std::to_string(gl) + "_sl" + + std::to_string(sl) + "_hpl_y" + mtag; + auto* env = makeEnv(n, 0.5 * hplZ, zCursor + 0.5 * hplZ); + CaloFibreHPLayer::build(env, alMat, psMat, n, 0.0, iHPL, cfg.plate_xy_mm, + cfg.hpl_thickness_mm, cfg.fiber_diameter_mm, + cfg.fiber_core_diameter_mm, true, mtag); + zCursor += hplZ; + ++iHPL; + ++gl; + ++sl; + + } else if (code == 6) { + const std::string n = "/SHiP/calorimeter/hcal/gl" + std::to_string(gl) + "_sl" + + std::to_string(sl) + "_hpl_x" + mtag; + auto* env = makeEnv(n, 0.5 * hplZ, zCursor + 0.5 * hplZ); + CaloFibreHPLayer::build(env, alMat, psMat, n, 0.0, iHPL, cfg.plate_xy_mm, + cfg.hpl_thickness_mm, cfg.fiber_diameter_mm, + cfg.fiber_core_diameter_mm, false, mtag); + zCursor += hplZ; + ++iHPL; + ++gl; + ++sl; + + } else if (code == 8) { + zCursor += airGapZ; + ++gl; + + } else { + throw std::runtime_error("CalorimeterFactory: unknown layer code in 'layers2': " + + std::to_string(code)); + } + } } } // namespace SHiPGeometry diff --git a/subsystems/Calorimeter/test_calorimeter.cpp b/subsystems/Calorimeter/test_calorimeter.cpp index 7ad72f0..456a304 100644 --- a/subsystems/Calorimeter/test_calorimeter.cpp +++ b/subsystems/Calorimeter/test_calorimeter.cpp @@ -1,6 +1,7 @@ // SPDX-License-Identifier: LGPL-3.0-or-later // Copyright (C) CERN for the benefit of the SHiP Collaboration +#include "Calorimeter/CalorimeterConfig.h" #include "Calorimeter/CalorimeterFactory.h" #include "SHiPGeometry/SHiPMaterials.h" @@ -9,18 +10,48 @@ #include #include +#include +using SHiPGeometry::CalorimeterConfig; +using SHiPGeometry::CalorimeterFactory; using SHiPGeometry::SHiPMaterials; TEST_CASE("CalorimeterBuilds", "[calorimeter]") { SHiPMaterials materials; - SHiPGeometry::CalorimeterFactory factory(materials); + CalorimeterFactory factory(materials); GeoPhysVol* calo = factory.build(); REQUIRE(calo != nullptr); auto* box = dynamic_cast(calo->getLogVol()->getShape()); REQUIRE(box != nullptr); - // Container: 3.00×3.50×1.45 m half-sizes + // Fixed envelope: 3.00 × 3.50 × 1.45 m half-sizes CHECK(box->getXHalfLength() == 3000.0); CHECK(box->getYHalfLength() == 3500.0); CHECK(box->getZHalfLength() == 1450.0); } + +TEST_CASE("CalorimeterHasChildren", "[calorimeter]") { + SHiPMaterials materials; + CalorimeterFactory factory(materials); + GeoPhysVol* calo = factory.build(); + REQUIRE(calo != nullptr); + // The container must have at least one child (ECAL layers + HCAL layers) + CHECK(calo->getNChildVols() >= 1u); +} + +TEST_CASE("TotalStackZPositive", "[calorimeter]") { + // Minimal config exercising all layer types + CalorimeterConfig cfg; + cfg.layers = {7, 1, 7, 2}; // lead, PVT-H, lead, PVT-V + cfg.layers2 = {7, 1}; // iron, PVT-H + cfg.lead_thickness_mm = 3.0; + cfg.scint_thickness_mm = 10.0; + cfg.iron_thickness_mm = 10.0; + cfg.hpl_thickness_mm = 10.0; + cfg.fiber_diameter_mm = 1.2; + cfg.fiber_core_diameter_mm = 1.0; + cfg.gap_ecal_hcal_mm = 0.0; + // layers: 2×lead + 2×scint = 6+20 = 26 mm + // layers2: 1×iron + 1×scint = 10+10 = 20 mm + // total = 46 mm + CHECK_THAT(CalorimeterFactory::totalStackZ(cfg), Catch::Matchers::WithinAbs(46.0, 1e-9)); +} diff --git a/subsystems/Trackers/README.md b/subsystems/Trackers/README.md index ce08c1f..50a89df 100644 --- a/subsystems/Trackers/README.md +++ b/subsystems/Trackers/README.md @@ -4,40 +4,101 @@ Straw tube tracking stations. ## Description -The Trackers subsystem implements 4 tracking stations for the SHiP spectrometer. Each station currently consists of an empty air box at the correct z-position. The full implementation requires straw tube modules with individual straws, stereo views, and support structures. +The Trackers subsystem implements 4 tracking stations for the SHiP spectrometer. +Stations 1–2 sit upstream of the main spectrometer magnet, stations 3–4 +downstream. Each station is a stack of 4 stereo views (layers); each view +contains a hollow aluminium frame and two staggered sub-layers of straw tubes. + +Geometry parameters (straw radius, length, pitch, stereo angle, frame size) +were ported from the standalone `StrawTrackerBuilder` reference; the per-station +envelope and station Z positions are taken from `subsystem_envelopes.csv`. ## Geometry Tree ``` -TrackersContainer (Air, 6000×6860×6000 mm) - ├─ TrackerStation_1 (Air, 6000×6860×1000 mm) z = 84070 mm - ├─ TrackerStation_2 (Air) z = 86070 mm - ├─ TrackerStation_3 (Air) z = 93070 mm - └─ TrackerStation_4 (Air) z = 95070 mm +/SHiP/trackers (Air, 6000 × 6860 × 12000 mm) + ├─ /SHiP/trackers/station_1 (Air, 6000 × 6860 × 1000 mm) z = 84070 mm + │ └─ Layer_j [j = 0..3] air box, rotated by ±2.3° about Z + │ ├─ StrawFrame hollow rectangle (Aluminium, GeoShapeSubtraction) + │ │ outer 4210 × 6230 mm, aperture 4010 × 6030 mm + │ ├─ StrawSubLayer_0 air slab at z = -10.55 mm (nominal) + │ │ └─ Straw [i = 0..299] wall (Mylar) + gas (Ar/CO₂) tubes + │ └─ StrawSubLayer_1 air slab at z = +10.55 mm (Y-staggered) + │ └─ Straw [i = 0..299] + ├─ /SHiP/trackers/station_2 z = 86070 mm + ├─ /SHiP/trackers/station_3 z = 93070 mm + └─ /SHiP/trackers/station_4 z = 95070 mm ``` -Position in world: centred at z = 89570 mm (average of stations 1 and 4). -Stations 1-2 are upstream of the magnet, stations 3-4 downstream. +Position in world: container centred at z = 89570 mm +(`(s_station1Z + s_station4Z) / 2`). + +### Straw specification + +| Parameter | Value | +|------------------|--------------------------------------| +| Outer radius | 10 mm (20 mm diameter) | +| Length | 4000 mm (along local X) | +| Wall thickness | 30 µm Mylar | +| Fill gas | Ar/CO₂ 70/30 by mass | +| Pitch | 20 mm | +| Straws/sub-layer | 300 | +| Sub-layers/view | 2 (second is +½-pitch staggered in Y)| +| Stereo angles | +2.3°, −2.3°, +2.3°, −2.3° | + +The straw axis is rotated +90° about Y so that the (locally Z-aligned) +`GeoTube` ends up pointing along the layer's X axis. ## Materials -| Material | Density | Usage | -|----------|-------------|----------------------| -| Air | 1.29 mg/cm³ | Container & stations | +| Material | Density | Composition | Usage | +|---------------|------------------|--------------------------------------------|----------------------| +| Air | 1.29 mg/cm³ | N 75.5%, O 23.1%, Ar 1.4% | Container, sub-layers | +| Aluminium | 2.70 g/cm³ | Al | View frames | +| Mylar | 1.39 g/cm³ | C₁₀H₈O₄ (PET, mass-fractioned) | Straw walls | +| ArCO2_70_30 | 1.56 × 10⁻³ g/cm³| Ar 70%, CO₂ 30% (by mass) | Straw gas (sensitive)| + +`Mylar` and `ArCO2_70_30` are added to `SHiPMaterials::createMaterials()` +in this PR; the elements they require (C, H, O, Ar) were already present. + +## Implementation notes + +### Shared `GeoLogVol`s + +The factory builds **one** `GeoLogVol` per repeated geometric element (straw +wall, straw gas, sub-layer envelope ×2, layer envelope, frame) and reuses it +across every placement. The 9600 straw placements therefore generate only two +straw `GeoLogVol`s instead of two per placement. + +### Stereo rotation + +Each layer envelope is built axis-aligned (straws along local X). The +`±2.3°` rotation about Z is applied at placement time inside +`buildStation()`, so the frame and both sub-layers rotate together with the +view — matching the physical reality where the frame holds the straws. + +### Solid-wall straws + +The straw wall is a **solid** Mylar `GeoTube` with the gas tube as its single +daughter. This avoids the mother/daughter overlap that hollow walls plus an +inner gas cylinder produce, which Geant4's overlap checker flags on every +straw. Physics-wise the result is identical: the gas material is used inside +the gas volume and Mylar everywhere else. ## Status -- [x] C++ implementation (envelope only) -- [ ] Implement straw tube geometry -- [ ] Add stereo views and support structures +- [x] C++ implementation (full straw geometry) +- [x] Stereo views and material frames +- [x] Mylar walls and Ar/CO₂ gas materials +- [ ] Anode wires (currently omitted — see TODO) - [ ] Verification against GDML +- [ ] Cross-check with FairShip's `strawtubesGeo` ## TODO -- Implement straw tube modules within each station (major work) - - Individual straw tubes (mylar + gas) - - 4 views per station (Y, U, V, Y') with stereo angles - - Support frames and service volumes -- Add straw tube gas material (Ar/CO2 mixture) to SHiPMaterials -- Add mylar material to SHiPMaterials -- Verify station positions against GDML reference +- Add a thin tungsten wire (~25 µm) on the gas tube axis if downstream + simulation needs an explicit anode. +- Verify station positions and layer pitch against the GDML reference once + it's available. +- Consider config-driving the geometry (analogous to `calo.cfg`) once a + second tracker variant lands. diff --git a/subsystems/Trackers/include/Trackers/TrackersFactory.h b/subsystems/Trackers/include/Trackers/TrackersFactory.h index 51551ea..7b7ba9d 100644 --- a/subsystems/Trackers/include/Trackers/TrackersFactory.h +++ b/subsystems/Trackers/include/Trackers/TrackersFactory.h @@ -3,50 +3,104 @@ #pragma once +#include + +class GeoLogVol; class GeoPhysVol; +class GeoVPhysVol; namespace SHiPGeometry { class SHiPMaterials; /** - * @brief Factory for the Trackers (straw tube tracking stations) geometry + * @brief Factory for the Trackers (straw tube tracking stations) geometry. + * + * Builds 4 stations of straw tubes spanning the SHiP spectrometer, each station + * containing four stereo views (layers) of double-staggered straws inside an + * aluminium support frame. + * + * Layout (per station, after stereo rotation about Z): + * - 4 layers : views with stereo angles +/-2.3 deg, alternating + * - 1 frame : hollow rectangle (outer - aperture, GeoShapeSubtraction) + * - 2 sub-layers (nominal + half-pitch staggered) + * - 300 straws each, 20 mm pitch in Y, straw axis along X (4 m long) + * - mylar wall 30 um, Ar/CO2 70/30 gas + * + * Station Z positions (centres, mm) match subsystem_envelopes.csv: + * station 1: 84070 (z = 83.57-84.57 m) + * station 2: 86070 + * station 3: 93070 (downstream of magnet) + * station 4: 95070 * - * Creates 4 tracking stations from GDML reference (statbox solid): - * - Station 1: Z 83.57-84.57 m → centre 84.07 m - * - Station 2: Z 85.57-86.57 m → centre 86.07 m - * - Station 3: Z 92.57-93.57 m → centre 93.07 m - * - Station 4: Z 94.57-95.57 m → centre 95.07 m - * GDML statbox: x=600 cm, y=686 cm, z=100 cm → half: 300×343×50 cm + * Container envelope (mm half-sizes) and per-station envelope are unchanged + * from the previous placeholder, so test_trackers' bound checks still apply. */ class TrackersFactory { public: explicit TrackersFactory(SHiPMaterials& materials); ~TrackersFactory() = default; - /** - * @brief Build the Trackers geometry - * @return Pointer to container volume with all 4 stations - */ + /** Build and return the tracker container volume with 4 populated stations. */ GeoPhysVol* build(); - private: - SHiPMaterials& m_materials; - - // Station dimensions from GDML statbox (mm) + // ── Per-station envelope (unchanged from the placeholder, GDML statbox) ── static constexpr double s_halfX = 3000.0; // 300 cm - static constexpr double s_halfY = 3430.0; // 343 cm (GDML y=686 cm) - static constexpr double s_halfZ = 500.0; // 50 cm + static constexpr double s_halfY = 3430.0; // 343 cm + static constexpr double s_halfZ = 500.0; // 50 cm - // Station Z positions (centres, in mm from origin) - static constexpr double s_station1Z = 84070.0; // 84.07 m - static constexpr double s_station2Z = 86070.0; // 86.07 m - static constexpr double s_station3Z = 93070.0; // 93.07 m - static constexpr double s_station4Z = 95070.0; // 95.07 m + // Station Z centres (mm from world origin). + static constexpr double s_station1Z = 84070.0; + static constexpr double s_station2Z = 86070.0; + static constexpr double s_station3Z = 93070.0; + static constexpr double s_station4Z = 95070.0; - // Container dimensions (spans all stations) + // Container spans station 1 -> station 4 (centre on average of extremes). static constexpr double s_containerHalfZ = (s_station4Z - s_station1Z) / 2.0 + s_halfZ; static constexpr double s_containerCentreZ = (s_station1Z + s_station4Z) / 2.0; + + // ── Straw geometry (mm, deg) ───────────────────────────────────────────── + static constexpr int s_nStations = 4; + static constexpr int s_nLayers = 4; // stereo views per station + static constexpr int s_nSubLayers = 2; // staggered pair per layer + static constexpr double s_strawRadius = 10.0; // 1 cm radius + static constexpr double s_strawLength = 4000.0; + static constexpr double s_wallThick = 0.030; // 30 um Mylar + static constexpr double s_apertureX = 4000.0; + static constexpr double s_apertureY = 6000.0; + static constexpr int s_nStraws = static_cast(s_apertureY / (2.0 * s_strawRadius)); // 300 + static constexpr double s_stereoAngleDeg = 2.3; + static constexpr double s_frameWidthX = 100.0; + static constexpr double s_frameWidthY = 100.0; + static constexpr double s_frameHalfZ = 22.0; + + private: + SHiPMaterials& m_materials; + + // Cached LogVols built once per build() and reused across all placements + // (one wall + one gas LogVol covers all 9600 straws, one shared LogVol per + // sub-layer variant, etc. — keeps the in-memory tree small). + GeoLogVol* m_strawWallLog = nullptr; + GeoLogVol* m_strawGasLog = nullptr; + GeoLogVol* m_subLayerNominal = nullptr; + GeoLogVol* m_subLayerShifted = nullptr; + GeoLogVol* m_layerLog = nullptr; + GeoLogVol* m_frameLog = nullptr; + + // ── Internal builders ──────────────────────────────────────────────────── + /** Build the LogVols that are shared across stations (straws, sub-layers, + * layer envelope, frame). Called once at the top of build(). */ + void buildSharedLogVols(); + + /** Build one station volume: 4 stereo-rotated layers in air. */ + GeoPhysVol* buildStation(int stationIndex); + + /** Place one stereo layer (frame + 2 sub-layers) inside @p station. */ + void placeLayer(GeoPhysVol* station, int layerIndex, double signedAngleDeg) const; + + /** Place one sub-layer of @p s_nStraws straws inside @p layer. + * @param shifted if true, the straws are staggered by +1/2 pitch in Y. */ + void placeSubLayer(GeoPhysVol* layer, bool shifted) const; }; } // namespace SHiPGeometry diff --git a/subsystems/Trackers/src/TrackersFactory.cpp b/subsystems/Trackers/src/TrackersFactory.cpp index 6171857..d001e1d 100644 --- a/subsystems/Trackers/src/TrackersFactory.cpp +++ b/subsystems/Trackers/src/TrackersFactory.cpp @@ -11,42 +11,252 @@ #include #include #include +#include #include +#include +#include +#include #include namespace SHiPGeometry { +namespace { + +// ── Internal layout helpers (private to this TU) ───────────────────────────── +// +// The straw pattern fills a (s_apertureX x s_apertureY) rectangle. Around it +// sits a hollow material frame; outside that, an air "layer envelope" that +// gets stereo-rotated by the parent station. All numbers are in millimetres. +// +// Aperture clearance: extra space around the straw pattern so that the +// sub-layer envelopes — which extend by +/- one straw radius in Z and (for +// the staggered sub-layer) +1/2 pitch in Y — fit inside the frame without +// touching its walls. +constexpr double kApClearX = 5.0; +constexpr double kApClearY = 15.0; + +// Small slack between the layer envelope and the frame outer surface (and +// between the sub-layer envelope and the frame aperture). Without this, +// CheckOverlaps in Geant4 flags cosmetic touches. +constexpr double kEnvClearance = 5.0; + +// Per-straw additional Z separation between the two sub-layer envelopes so +// they do not share a face at z = 0 in the layer frame. +constexpr double kSubLayerZSlack = 0.55; + +// Frame aperture (inner hole) half-sizes in mm. +constexpr double kApHalfX = TrackersFactory::s_apertureX / 2.0 + kApClearX; // 2005 +constexpr double kApHalfY = TrackersFactory::s_apertureY / 2.0 + kApClearY; // 3015 + +// Frame outer half-sizes in mm. +constexpr double kFrHalfX = kApHalfX + TrackersFactory::s_frameWidthX; // 2105 +constexpr double kFrHalfY = kApHalfY + TrackersFactory::s_frameWidthY; // 3115 + +// Layer envelope half-sizes in mm. +constexpr double kLayHalfX = kFrHalfX + kEnvClearance; // 2110 +constexpr double kLayHalfY = kFrHalfY + kEnvClearance; // 3120 +constexpr double kLayHalfZ = TrackersFactory::s_frameHalfZ + kEnvClearance; // 27 + +// Sub-layer envelope half-sizes in mm. Slightly smaller than the aperture so +// it sits cleanly inside the frame; thick enough in Z to wrap one straw. +constexpr double kFrameClearance = 0.5; +constexpr double kSlHalfX = kApHalfX - kFrameClearance; // 2004.5 +constexpr double kSlHalfY = kApHalfY - kFrameClearance; // 3014.5 +constexpr double kSlHalfZ = TrackersFactory::s_strawRadius + 0.5; // 10.5 + +// Z stack of layers within a station: kNLayers air slabs of half-thickness +// kLayHalfZ separated by a small gap. +constexpr double kLayerGap = 5.0; +constexpr double kLayerPitch = 2.0 * kLayHalfZ + kLayerGap; // 59 + +double signedStereoDeg(int layerIndex) { + // Layers 0,2 -> +angle (u view); layers 1,3 -> -angle (v view). + return (layerIndex % 2 == 0 ? +1.0 : -1.0) * TrackersFactory::s_stereoAngleDeg; +} + +double layerZInStation(int layerIndex) { + // Centre the 4-layer stack on z = 0 within the station envelope. + return -0.5 * (TrackersFactory::s_nLayers - 1) * kLayerPitch + layerIndex * kLayerPitch; +} + +} // namespace + TrackersFactory::TrackersFactory(SHiPMaterials& materials) : m_materials(materials) {} +// ───────────────────────────────────────────────────────────────────────────── +// build() +// ───────────────────────────────────────────────────────────────────────────── GeoPhysVol* TrackersFactory::build() { const GeoMaterial* air = m_materials.requireMaterial("Air"); - // Create container volume that spans all 4 stations + buildSharedLogVols(); + + // Tracker container that spans all 4 stations. Half-Z is large enough to + // include the most upstream and downstream station envelopes. auto* containerBox = new GeoBox(s_halfX, s_halfY, s_containerHalfZ); auto* containerLog = new GeoLogVol("/SHiP/trackers", containerBox, air); auto* containerPhys = new GeoPhysVol(containerLog); - // Create and place individual stations - const double stationZ[4] = {s_station1Z, s_station2Z, s_station3Z, s_station4Z}; + constexpr std::array stationZ = {s_station1Z, s_station2Z, s_station3Z, s_station4Z}; - for (int i = 0; i < 4; ++i) { - auto* stationBox = new GeoBox(s_halfX, s_halfY, s_halfZ); - std::string stationName = "/SHiP/trackers/station_" + std::to_string(i + 1); - auto* stationLog = new GeoLogVol(stationName, stationBox, air); - auto* stationPhys = new GeoPhysVol(stationLog); + for (int i = 0; i < s_nStations; ++i) { + GeoPhysVol* stationPhys = buildStation(i); - // Position relative to container centre - double relativeZ = stationZ[i] - s_containerCentreZ; - GeoTrf::Transform3D stationTrf = GeoTrf::Translate3D(0.0, 0.0, relativeZ); + const std::string stationName = "/SHiP/trackers/station_" + std::to_string(i + 1); + const double relativeZ = stationZ[i] - s_containerCentreZ; containerPhys->add(new GeoNameTag(stationName)); - containerPhys->add(new GeoIdentifierTag(i)); - containerPhys->add(new GeoTransform(stationTrf)); + containerPhys->add(new GeoIdentifierTag(i + 1)); + containerPhys->add(new GeoTransform(GeoTrf::Translate3D(0.0, 0.0, relativeZ))); containerPhys->add(stationPhys); } return containerPhys; } +// ───────────────────────────────────────────────────────────────────────────── +// buildSharedLogVols() +// +// Builds every LogVol that is repeated across stations exactly once, so the +// final tree contains O(few) LogVols instead of O(9600). Placement-time +// individuation is done with GeoNameTag and GeoIdentifierTag. +// ───────────────────────────────────────────────────────────────────────────── +void TrackersFactory::buildSharedLogVols() { + const GeoMaterial* air = m_materials.requireMaterial("Air"); + const GeoMaterial* mylar = m_materials.requireMaterial("Mylar"); + const GeoMaterial* gasArCO2 = m_materials.requireMaterial("ArCO2_70_30"); + const GeoMaterial* alu = m_materials.requireMaterial("Aluminium"); + + // ── Straw wall + gas (one pair shared by every straw in every layer) ──── + // The wall is built as a SOLID Mylar tube; the gas tube sits inside it as + // a daughter and physics-wise replaces the wall material in its volume. + // This avoids the mother/daughter overlap that the (rMin, rMax) hollow + // tube + gas-cylinder combination triggers in Geant4's overlap checker. + const double rGas = s_strawRadius - s_wallThick; + const double rWall = s_strawRadius; + const double half = 0.5 * s_strawLength; + + auto* wallTube = new GeoTube(0.0, rWall, half); + m_strawWallLog = new GeoLogVol("StrawWall", wallTube, mylar); + + auto* gasTube = new GeoTube(0.0, rGas, half); + m_strawGasLog = new GeoLogVol("StrawGas", gasTube, gasArCO2); + + // ── Sub-layer envelopes (nominal + shifted) ───────────────────────────── + auto* slBox = new GeoBox(kSlHalfX, kSlHalfY, kSlHalfZ); + m_subLayerNominal = new GeoLogVol("StrawSubLayer_nominal", slBox, air); + m_subLayerShifted = new GeoLogVol("StrawSubLayer_shifted", slBox, air); + + // ── Layer envelope (one shape, reused for all 16 placements) ──────────── + auto* layBox = new GeoBox(kLayHalfX, kLayHalfY, kLayHalfZ); + m_layerLog = new GeoLogVol("StrawLayer", layBox, air); + + // ── Material frame: outer rectangle minus inner aperture ──────────────── + auto* outerBox = new GeoBox(kFrHalfX, kFrHalfY, s_frameHalfZ); + auto* innerBox = new GeoBox(kApHalfX, kApHalfY, s_frameHalfZ + 1.0); + auto* frameShape = new GeoShapeSubtraction(outerBox, innerBox); + m_frameLog = new GeoLogVol("StrawFrame", frameShape, alu); +} + +// ───────────────────────────────────────────────────────────────────────────── +// buildStation() +// ───────────────────────────────────────────────────────────────────────────── +GeoPhysVol* TrackersFactory::buildStation(int stationIndex) { + const GeoMaterial* air = m_materials.requireMaterial("Air"); + + // Each station uses the per-station envelope from subsystem_envelopes.csv. + // The internal stack (4 layers, ~125 mm half-Z) sits comfortably inside + // the 500 mm half-Z envelope — leaving room for service material later. + // The station LogVol name is unique per station so test_trackers' lookup + // by /SHiP/trackers/station_ continues to work. + const std::string stationName = "/SHiP/trackers/station_" + std::to_string(stationIndex + 1); + + auto* stationBox = new GeoBox(s_halfX, s_halfY, s_halfZ); + auto* stationLog = new GeoLogVol(stationName, stationBox, air); + auto* stationPhys = new GeoPhysVol(stationLog); + + for (int j = 0; j < s_nLayers; ++j) { + placeLayer(stationPhys, j, signedStereoDeg(j)); + } + return stationPhys; +} + +// ───────────────────────────────────────────────────────────────────────────── +// placeLayer() +// ───────────────────────────────────────────────────────────────────────────── +void TrackersFactory::placeLayer(GeoPhysVol* station, int layerIndex, double signedAngleDeg) const { + auto* layerPhys = new GeoPhysVol(m_layerLog); + + // ── 1. Frame at z = 0 (rotates with the view) ─────────────────────────── + layerPhys->add(new GeoNameTag("StrawFrame")); + layerPhys->add(new GeoIdentifierTag(100 + layerIndex)); + layerPhys->add(new GeoTransform(GeoTrf::Transform3D::Identity())); + layerPhys->add(new GeoPhysVol(m_frameLog)); + + // ── 2. Two sub-layers, staggered in Y and offset in Z ─────────────────── + const double dz = s_strawRadius + kSubLayerZSlack; // ~10.55 mm + + layerPhys->add(new GeoNameTag("StrawSubLayer_0")); + layerPhys->add(new GeoIdentifierTag(0)); + layerPhys->add(new GeoTransform(GeoTrf::TranslateZ3D(-dz))); + { + auto* sl0 = new GeoPhysVol(m_subLayerNominal); + placeSubLayer(sl0, /*shifted=*/false); + layerPhys->add(sl0); + } + + layerPhys->add(new GeoNameTag("StrawSubLayer_1")); + layerPhys->add(new GeoIdentifierTag(1)); + layerPhys->add(new GeoTransform(GeoTrf::TranslateZ3D(+dz))); + { + auto* sl1 = new GeoPhysVol(m_subLayerShifted); + placeSubLayer(sl1, /*shifted=*/true); + layerPhys->add(sl1); + } + + // ── 3. Place the (now populated) layer envelope inside the station ───── + const double angleRad = signedAngleDeg * std::numbers::pi / 180.0; + const double zLay = layerZInStation(layerIndex); + + GeoTrf::Transform3D xfLayer = GeoTrf::TranslateZ3D(zLay) * GeoTrf::RotateZ3D(angleRad); + + station->add(new GeoNameTag("Layer_" + std::to_string(layerIndex))); + station->add(new GeoIdentifierTag(layerIndex)); + station->add(new GeoTransform(xfLayer)); + station->add(layerPhys); +} + +// ───────────────────────────────────────────────────────────────────────────── +// placeSubLayer() +// +// Lays down s_nStraws straws along Y, pitch = 2 * straw radius. Each straw is +// a Mylar wall PhysVol with a gas PhysVol daughter; both LogVols are shared +// across every straw in the geometry. The straw axis is along the layer's +// local X (achieved by a +90 deg rotation about Y from the GeoTube's native +// Z axis). +// ───────────────────────────────────────────────────────────────────────────── +void TrackersFactory::placeSubLayer(GeoPhysVol* layer, bool shifted) const { + const double pitch = 2.0 * s_strawRadius; + const double yStart = -(s_nStraws - 1) * 0.5 * pitch; + const double yShift = shifted ? s_strawRadius : 0.0; + + for (int i = 0; i < s_nStraws; ++i) { + const double yStraw = yStart + i * pitch + yShift; + + // The wall PhysVol carries the gas as its single daughter. Constructing + // it here per straw is necessary because each placement gets its own + // GeoPhysVol node — but the underlying LogVols are shared. + auto* wallPhys = new GeoPhysVol(m_strawWallLog); + wallPhys->add(new GeoNameTag("StrawGas")); + wallPhys->add(new GeoPhysVol(m_strawGasLog)); + + layer->add(new GeoNameTag("Straw")); + layer->add(new GeoIdentifierTag(i)); + layer->add(new GeoTransform(GeoTrf::TranslateY3D(yStraw) * + GeoTrf::RotateY3D(std::numbers::pi / 2.0))); + layer->add(wallPhys); + } +} + } // namespace SHiPGeometry diff --git a/subsystems/Trackers/test_trackers.cpp b/subsystems/Trackers/test_trackers.cpp index f0d1c11..269ff5a 100644 --- a/subsystems/Trackers/test_trackers.cpp +++ b/subsystems/Trackers/test_trackers.cpp @@ -13,6 +13,7 @@ #include using SHiPGeometry::SHiPMaterials; +using SHiPGeometry::TrackersFactory; static const GeoVPhysVol* findChild(const GeoVPhysVol* parent, const std::string& name) { for (unsigned int i = 0; i < parent->getNChildVols(); ++i) { @@ -24,10 +25,11 @@ static const GeoVPhysVol* findChild(const GeoVPhysVol* parent, const std::string return nullptr; } -// CSV limits: Trackers per-station halfX ≤ 3000, halfY ≤ 3500, halfZ ≤ 500 +// CSV limits: Trackers per-station halfX <= 3000, halfY <= 3500, halfZ <= 500. +// Pre-existing test — kept verbatim so the new straw geometry doesn't regress it. TEST_CASE("TrackersWithinEnvelope", "[trackers]") { SHiPMaterials materials; - SHiPGeometry::TrackersFactory factory(materials); + TrackersFactory factory(materials); GeoPhysVol* tc = factory.build(); REQUIRE(tc != nullptr); const GeoVPhysVol* st1 = findChild(tc, "/SHiP/trackers/station_1"); @@ -39,3 +41,26 @@ TEST_CASE("TrackersWithinEnvelope", "[trackers]") { CHECK(box->getYHalfLength() <= 3500.0); CHECK(box->getZHalfLength() <= 500.0); } + +TEST_CASE("TrackersHaveAllFourStations", "[trackers]") { + SHiPMaterials materials; + TrackersFactory factory(materials); + GeoPhysVol* tc = factory.build(); + REQUIRE(tc != nullptr); + for (int i = 1; i <= 4; ++i) { + const std::string name = "/SHiP/trackers/station_" + std::to_string(i); + INFO("Missing station: " << name); + REQUIRE(findChild(tc, name) != nullptr); + } +} + +TEST_CASE("TrackerStationContainsLayers", "[trackers]") { + SHiPMaterials materials; + TrackersFactory factory(materials); + GeoPhysVol* tc = factory.build(); + REQUIRE(tc != nullptr); + const GeoVPhysVol* st = findChild(tc, "/SHiP/trackers/station_1"); + REQUIRE(st != nullptr); + // Each station now holds 4 stereo layers; previously this was empty. + CHECK(st->getNChildVols() == 4u); +} diff --git a/tests/test_consistency.cpp b/tests/test_consistency.cpp index ff1d270..9bff5b3 100644 --- a/tests/test_consistency.cpp +++ b/tests/test_consistency.cpp @@ -76,9 +76,9 @@ TEST_CASE("ConsistencyTest.ExpectedSubsystemCount", "[consistency]") { REQUIRE(world != nullptr); auto subsystems = collectSubsystems(world); - // 8 subsystems: target, muon_shield, upstream_tagger, decay_volume, + // 9 subsystems: target, muon_shield, upstream_tagger, decay_volume, sbt, // trackers, magnet, timing_detector, calorimeter - CHECK(subsystems.size() == 8u); // NOLINT(readability/check) + CHECK(subsystems.size() == 9u); // NOLINT(readability/check) } TEST_CASE("ConsistencyTest.SubsystemsGenerallyInZOrder", "[consistency]") { @@ -110,9 +110,12 @@ TEST_CASE("ConsistencyTest.NoUnexpectedZOverlaps", "[consistency]") { // The trackers container intentionally spans across the magnet // (stations 1-2 before, stations 3-4 after), so that pair is allowed to overlap. + // The SBT wraps the decay volume (same Z range by design), so they also overlap. auto isAllowedOverlap = [](const std::string& a, const std::string& b) { return (a == "/SHiP/trackers" && b == "/SHiP/magnet") || - (a == "/SHiP/magnet" && b == "/SHiP/trackers"); + (a == "/SHiP/magnet" && b == "/SHiP/trackers") || + (a == "/SHiP/sbt" && b == "/SHiP/decay_volume") || + (a == "/SHiP/decay_volume" && b == "/SHiP/sbt"); }; // Sort by Z centre