From 05723216b45c094ffb34ef5724e513b50e55d842 Mon Sep 17 00:00:00 2001 From: Nicola Nicassio Date: Fri, 17 Jul 2026 19:57:32 +0200 Subject: [PATCH 1/4] Latest v3b.1 RICH geometry with quadrants, modules and shielding --- .../base/include/RICHBase/RICHBaseParam.h | 93 ++- .../include/RICHSimulation/RICHRing.h | 100 +-- .../ALICE3/RICH/simulation/src/Detector.cxx | 708 +++++++++++++++++- .../ALICE3/RICH/simulation/src/RICHRing.cxx | 653 +++++++++++++--- 4 files changed, 1331 insertions(+), 223 deletions(-) diff --git a/Detectors/Upgrades/ALICE3/RICH/base/include/RICHBase/RICHBaseParam.h b/Detectors/Upgrades/ALICE3/RICH/base/include/RICHBase/RICHBaseParam.h index a9f2f7fbba5d1..52f5725ed032a 100644 --- a/Detectors/Upgrades/ALICE3/RICH/base/include/RICHBase/RICHBaseParam.h +++ b/Detectors/Upgrades/ALICE3/RICH/base/include/RICHBase/RICHBaseParam.h @@ -20,34 +20,89 @@ namespace o2 namespace rich { struct RICHBaseParam : public o2::conf::ConfigurableParamHelper { - float zBaseSize = 18.6; // cm (18.4 in v3) - float rMax = 131.0; // cm (117.0 in v3) - float rMin = 104.0; // cm (90.0 in v3) - float radiatorThickness = 2.0; // cm - float detectorThickness = 0.2; // cm - float zRichLength = 700.0; // cm + double zBaseSize = 18.6; // cm (18.4 in v3) + double rMax = 131.0; // cm (117.0 in v3) + double rMin = 104.0; // cm (90.0 in v3) + double radiatorThickness = 2.0; // cm + double zRichLength = 700.0; // cm int nRings = 11; // (25 in v3) int nTiles = 44; // (36 in v3) bool oddGeom = true; // (false in v3) + + // The active and passive silicon thicknesses must sum to detectorThickness. + double siliconeLayerThickness = 0.010; // cm: 0.1 mm resin layer in front + double detectorThickness = 0.1; // cm + double activeSiliconThickness = 0.01; // cm: 0.1 mm sensitive silicon + //double passiveSiliconThickness = 0.09f; // cm: (detectorThickness - activeSiliconThickness) + + //cylindrical aerogel layout + bool useCylindricalAerogel = true; + double cylindricalAerogelEtaRef = 0.85; + + // Enable geometry with rectangular modules + bool useRectangularModules = true; + + // Barrel photosensor active area. + double sipmActiveSizeZ = 18.0; // cm + double sipmActiveSizeRPhi = 17.0; // cm - // FWD and BWD RICH - bool enableFWDRich = false; - bool enableBWDRich = false; + // Gas refractive index (then scaled with chromaticity) + double nGasEffective = 1.0006; - float rFWDMin = 13.7413f; - float rFWDMax = 103.947f; + // Aerogel refractive index (then scaled with chromaticity) + double nAerogelEffective = 1.03; - // Aerogel: - float zAerogelMin = 375.f; - float zAerogelMax = 377.f; + // Parameters for geometry with quadrants + bool flagUseQuadrants = false; + // Opening between adjacent vessel quadrants, measured as a chord at shieldRMin. + double vesselPhiGap = 1.0; // cm + // Thickness of each lateral insulating wall at a quadrant boundary. + double vesselThicknessShieldingLateral = 1.0; // cm + // Rectangular size could be smaller with quadrants (< 17 cm depending on wall thickness) + double quadrantModuleSizeRPhi = 16.5; // cm - // Argon: - float zArgonMin = 377.f; - float zArgonMax = 407.f; + // Readout stack behind each SiPM plane, thicknesses along the local outward normal. + double pcb1Thickness = 0.4; // cm + double coolingPlateThickness = 0.4; // cm + double pcb2Thickness = 0.4; // cm + double pcb3Thickness = 0.4; // cm + // Surface-to-surface gaps between consecutive layers. + double gapSiPMToPCB1 = 0.10; // cm + double gapPCB1ToCoolingPlate = 0.10; // cm + double gapCoolingPlateToPCB2 = 0.10; // cm + double gapPCB2ToPCB3 = 0.10; // cm + + // Minimum edge-to-edge clearances used to avoid exact contacts between adjacent modules. + double moduleClearanceZ = 0.02; // cm + double moduleClearanceRPhi = 0.02; // cm + + // Shielding: + // Radial boundaries of the complete cylindrical enclosure. + double shieldRMin = 100.0; + double shieldRMax = 136.0; + // Radial thickness of the inner insulating wall. + double innerWallThickness = 2.0; + // Radial thickness of the outer insulating wall. + double outerWallThickness = 2.0; + // Full longitudinal length of the cylindrical side walls. + double shieldLengthZ = 220.0; + // Thickness of each insulating end cap along Z. + double endCapThicknessZ = 2.0; + // FWD and BWD RICH (legacy) + bool enableFWDRich = false; + bool enableBWDRich = false; + double rFWDMin = 13.7413; + double rFWDMax = 103.947; + // Aerogel: + double zAerogelMin = 375.; + double zAerogelMax = 377.; + // Argon: + double zArgonMin = 377.; + double zArgonMax = 407.; // Detector: - float zSiliconMin = 407.f; - float zSiliconMax = 407.2f; + double zSiliconMin = 407.; + double zSiliconMax = 407.2; O2ParamDef(RICHBaseParam, "RICHBase"); }; diff --git a/Detectors/Upgrades/ALICE3/RICH/simulation/include/RICHSimulation/RICHRing.h b/Detectors/Upgrades/ALICE3/RICH/simulation/include/RICHSimulation/RICHRing.h index 296e24cbd8f06..18c4c4e2f8d71 100644 --- a/Detectors/Upgrades/ALICE3/RICH/simulation/include/RICHSimulation/RICHRing.h +++ b/Detectors/Upgrades/ALICE3/RICH/simulation/include/RICHSimulation/RICHRing.h @@ -35,20 +35,20 @@ class Ring // z_ph: z position of the photosensitive surface (from the center) Ring(int rPosId, int nTilesPhi, - float rMin, - float rMax, - float radThick, - float radYmin, - float radYmax, - float radZ, - float photThick, - float photYmin, - float photYmax, - float photZ, - float radRad0, - float photRad0, - float aerDetDistance, - float thetaB, + double rMin, + double rMax, + double radThick, + double radYmin, + double radYmax, + double radZ, + double photThick, + double photYmin, + double photYmax, + double photZ, + double radRad0, + double photRad0, + double aerDetDistance, + double thetaB, const std::string motherName = "RICHV"); ~Ring() = default; @@ -60,10 +60,10 @@ class Ring private: int mPosId; // id of the ring int mNTiles; // number of modules - float mRRad; // max distance for radiators - float mRPhot; // max distance for photosensitive surfaces - float mRadThickness; // thickness of the radiator - float mPhotThickness; // thickness of the photosensitive surface + double mRRad; // max distance for radiators + double mRPhot; // max distance for photosensitive surfaces + double mRadThickness; // thickness of the radiator + double mPhotThickness; // thickness of the photosensitive surface ClassDef(Ring, 0); }; @@ -74,32 +74,32 @@ class FWDRich public: FWDRich() = default; FWDRich(std::string name, - float rMin, - float rMax, - float zAerogelMin, - float dZAerogel, - float zArgonMin, - float dZArgon, - float zSiliconMin, - float dZSilicon); + double rMin, + double rMax, + double zAerogelMin, + double dZAerogel, + double zArgonMin, + double dZArgon, + double zSiliconMin, + double dZSilicon); void createFWDRich(TGeoVolume* motherVolume); protected: std::string mName; - float mRmin; - float mRmax; + double mRmin; + double mRmax; // Aerogel: - float mZAerogelMin; - float mDZAerogel; + double mZAerogelMin; + double mDZAerogel; // Argon: - float mZArgonMin; - float mDZArgon; + double mZArgonMin; + double mDZArgon; // Silicon: - float mZSiliconMin; - float mDZSilicon; + double mZSiliconMin; + double mDZSilicon; ClassDef(FWDRich, 0); }; @@ -109,32 +109,32 @@ class BWDRich public: BWDRich() = default; BWDRich(std::string name, - float rMin, - float rMax, - float zAerogelMin, - float dZAerogel, - float zArgonMin, - float dZArgon, - float zSiliconMin, - float dZSilicon); + double rMin, + double rMax, + double zAerogelMin, + double dZAerogel, + double zArgonMin, + double dZArgon, + double zSiliconMin, + double dZSilicon); void createBWDRich(TGeoVolume* motherVolume); protected: std::string mName; - float mRmin; - float mRmax; + double mRmin; + double mRmax; // Aerogel: - float mZAerogelMin; - float mDZAerogel; + double mZAerogelMin; + double mDZAerogel; // Argon: - float mZArgonMin; - float mDZArgon; + double mZArgonMin; + double mDZArgon; // Silicon: - float mZSiliconMin; - float mDZSilicon; + double mZSiliconMin; + double mDZSilicon; ClassDef(BWDRich, 0); }; diff --git a/Detectors/Upgrades/ALICE3/RICH/simulation/src/Detector.cxx b/Detectors/Upgrades/ALICE3/RICH/simulation/src/Detector.cxx index 02719d6f93a00..fd415b76564de 100644 --- a/Detectors/Upgrades/ALICE3/RICH/simulation/src/Detector.cxx +++ b/Detectors/Upgrades/ALICE3/RICH/simulation/src/Detector.cxx @@ -15,6 +15,9 @@ #include #include #include +#include +#include +#include #include "DetectorsBase/Stack.h" #include "ITSMFTSimulation/Hit.h" @@ -27,6 +30,49 @@ namespace o2 { namespace rich { +namespace // quadrant equation solver +{ +double quadrantDeltaPhiEquation(double x, int nTilesPhi, double rMin, double totalBoundaryWidth) +{ + const double argument = totalBoundaryWidth * TMath::Cos(x / 2.0) / (2.0 * rMin); + if (TMath::Abs(argument) >= 1.0) { + return std::numeric_limits::quiet_NaN(); + } + const double rhs = 2.0 * TMath::Pi() / static_cast(nTilesPhi) - (8.0 / static_cast(nTilesPhi)) * TMath::ASin(argument); + return rhs - x; +} + +double solveQuadrantDeltaPhi(int nTilesPhi, double rMin, double totalBoundaryWidth) +{ + double lower = 0.0; + double upper = 1.1 * 2.0 * TMath::Pi() / static_cast(nTilesPhi); + double fLower = quadrantDeltaPhiEquation(lower, nTilesPhi, rMin, totalBoundaryWidth); + double fUpper = quadrantDeltaPhiEquation(upper, nTilesPhi, rMin, totalBoundaryWidth); + if (!std::isfinite(fLower) || !std::isfinite(fUpper) || fLower * fUpper > 0.0) { + return -1.0; + } + constexpr double tolerance = 1.0e-12; + constexpr int maxIterations = 200; + for (int iteration = 0; iteration < maxIterations; iteration++) { + const double middle = 0.5 * (lower + upper); + const double fMiddle = quadrantDeltaPhiEquation(middle, nTilesPhi, rMin, totalBoundaryWidth); + if (!std::isfinite(fMiddle)) { + return -1.0; + } + if (TMath::Abs(fMiddle) < tolerance || 0.5 * (upper - lower) < tolerance) { + return middle; + } + if (fLower * fMiddle < 0.0) { + upper = middle; + fUpper = fMiddle; + } else { + lower = middle; + fLower = fMiddle; + } + } + return 0.5 * (lower + upper); +} +} // namespace Detector::Detector() : o2::base::DetImpl("RCH", true), @@ -61,6 +107,10 @@ void Detector::ConstructGeometry() void Detector::createMaterials() { + auto& richPars = RICHBaseParam::Instance(); + const double nGasEffective = richPars.nGasEffective; + const double nAerogelEffective = richPars.nAerogelEffective; + int ifield = 2; // ? float fieldm = 10.0; // ? o2::base::Detector::initFieldTrackingParams(ifield, fieldm); @@ -89,12 +139,75 @@ void Detector::createMaterials() float epsilAerogel = 1.0E-4; // .10000E+01; float stminAerogel = 0.0; // cm "Default value used" + float tmaxfdCO2 = 0.1; // .10000E+01; // Degree + float stemaxCO2 = .10000E+01; // cm + float deemaxCO2 = 0.1; // 0.30000E-02; // Fraction of particle's energy 0SetCerenkov(globalMediumID(2, "AEROGEL"), nAerogelRindex, aerogelRindexEnergyGeV, aerogelAbsorptionOnRindexGrid, aerogelDetectionEfficiency, aerogelRindex); + // + // constexpr int nAerogelAbsorption = 2; + // double aerogelAbsorptionEnergyGeV[nAerogelAbsorption] = {1.0 * eVInGeV, 8.26561 * eVInGeV}; + // double aerogelAbsorptionLengthCm[nAerogelAbsorption] = {aerogelAbsorptionLengthCm, aerogelAbsorptionLengthCm}; + // mc->SetMaterialProperty(globalMediumID(2, "AEROGEL"), "ABSLENGTH", nAerogelAbsorption, aerogelAbsorptionEnergyGeV, aerogelAbsorptionLengthCm); + // + constexpr int nAerogelRayleigh = 22; + double aerogelRayleighEnergyGeV[nAerogelRayleigh] = {1.00 * eVInGeV, 1.06 * eVInGeV, 1.12 * eVInGeV, 1.18 * eVInGeV, 1.23984 * eVInGeV, 1.3051 * eVInGeV, 1.3776 * eVInGeV, 1.45864 * eVInGeV, 1.5498 * eVInGeV, 1.65312 * eVInGeV, 1.7712 * eVInGeV, 1.90745 * eVInGeV, 2.0664 * eVInGeV, 2.25426 * eVInGeV, 2.47968 * eVInGeV, 2.7552 * eVInGeV, 3.0996 * eVInGeV, 3.54241 * eVInGeV, 4.13281 * eVInGeV, 4.95937 * eVInGeV, 6.19921 * eVInGeV, 8.26561 * eVInGeV}; + double aerogelRayleighLengthCm[nAerogelRayleigh] = {543.253684, 430.307801, 345.247537, 280.204207, 229.885, 187.243, 150.828, 120.001, 94.1609, 72.7371, 55.1954, 41.0359, 29.7931, 21.0359, 14.3678, 9.42672, 5.88506, 3.44971, 1.86207, 0.897989, 0.367816, 0.116379}; + mc->SetMaterialProperty(globalMediumID(2, "AEROGEL"), "RAYLEIGH", nAerogelRayleigh, aerogelRayleighEnergyGeV, aerogelRayleighLengthCm); + + /// GAS + constexpr int nCO2Optical = 2; + double co2EnergyGeV[nCO2Optical] = {1.0 * eVInGeV, 8.26561 * eVInGeV}; + double co2Rindex[nCO2Optical] = {nGasEffective, nGasEffective}; // <- Target gas index for dielectrons + double co2AbsorptionLengthCm[nCO2Optical] = {1.0e5, 1.0e5}; + double co2DetectionEfficiency[nCO2Optical] = {0.0, 0.0}; + mc->SetCerenkov(globalMediumID(5, "CO2"), nCO2Optical, co2EnergyGeV, co2AbsorptionLengthCm, co2DetectionEfficiency, co2Rindex); + + /// SiO2 + constexpr int nSiO2Optical = 2; + double sio2EnergyGeV[nSiO2Optical] = {1.0 * eVInGeV, 8.26561 * eVInGeV}; + double sio2Rindex[nSiO2Optical] = {1.47, 1.47}; + double sio2AbsorptionLengthCm[nSiO2Optical] = {1.0e5, 1.0e5}; + double sio2DetectionEfficiency[nSiO2Optical] = {0.0, 0.0}; + mc->SetCerenkov(globalMediumID(9, "SIO2"), nSiO2Optical, sio2EnergyGeV, sio2AbsorptionLengthCm, sio2DetectionEfficiency, sio2Rindex); + + /// Silicone resin + constexpr int nSiliconeOptical = 2; + double siliconeEnergyGeV[nSiliconeOptical] = {1.0 * eVInGeV, 8.26561 * eVInGeV}; + double siliconeRindex[nSiliconeOptical] = {1.41, 1.41}; + double siliconeAbsorptionLengthCm[nSiliconeOptical] = {1.0e5, 1.0e5}; + double siliconeDetectionEfficiency[nSiliconeOptical] = {0.0, 0.0}; + mc->SetCerenkov(globalMediumID(10, "SILICONE"), nSiliconeOptical, siliconeEnergyGeV, siliconeAbsorptionLengthCm, siliconeDetectionEfficiency, siliconeRindex); + + /// Si (assuming same index as silicone resin as reflection losses are already included in PDE) + constexpr int nSiliconOptical = 2; + double siliconEnergyGeV[nSiliconOptical] = {1.0 * eVInGeV, 8.26561 * eVInGeV}; + double siliconRindex[nSiliconOptical] = {1.41, 1.41}; + double siliconAbsorptionLengthCm[nSiliconOptical] = {1.0e5, 1.0e5}; + double siliconDetectionEfficiency[nSiliconOptical] = {0.0, 0.0}; + mc->SetCerenkov(globalMediumID(3, "SILICON"), nSiliconOptical, siliconEnergyGeV, siliconAbsorptionLengthCm, siliconDetectionEfficiency, siliconRindex); + + // Si: outer layer just for photon absorption + constexpr int nSiliconAbsorberOptical = 2; + double siliconAbsorberEnergyGeV[nSiliconAbsorberOptical] = {1.0 * eVInGeV, 8.26561 * eVInGeV}; + double siliconAbsorberAbsorptionLengthCm[nSiliconAbsorberOptical] = {1.0e-7, 1.0e-7}; // 1 nm + mc->SetMaterialProperty(globalMediumID(11, "SILICON_ABSORBER"), "ABSLENGTH", nSiliconAbsorberOptical, siliconAbsorberEnergyGeV, siliconAbsorberAbsorptionLengthCm); } void Detector::createGeometry() @@ -145,8 +429,237 @@ void Detector::createGeometry() vRICH->SetTitle(vstrng); auto& richPars = RICHBaseParam::Instance(); + // Quadrant parameters + const bool flagUseQuadrants = richPars.flagUseQuadrants; + const double vesselPhiGap = richPars.vesselPhiGap; + const double vesselThicknessShieldingLateral = richPars.vesselThicknessShieldingLateral; + + // shielding parameters + double shieldRMin = richPars.shieldRMin; + double shieldRMax = richPars.shieldRMax; + double innerWallThickness = richPars.innerWallThickness; + double outerWallThickness = richPars.outerWallThickness; + double shieldLengthZ = richPars.shieldLengthZ; + double endCapThicknessZ = richPars.endCapThicknessZ; + + if (innerWallThickness <= 0.0 || outerWallThickness <= 0.0 || endCapThicknessZ <= 0.0 || shieldLengthZ <= 0.0) { + LOGP(fatal, "RICH shielding dimensions must be positive"); + } + + if (shieldRMin + innerWallThickness >= shieldRMax - outerWallThickness) { + LOGP(fatal, + "RICH shielding walls overlap: inner outer radius = {}, outer inner radius = {}", + shieldRMin + innerWallThickness, + shieldRMax - outerWallThickness); + } + + if (flagUseQuadrants) { + if (richPars.nTiles <= 0 || richPars.nTiles % 4 != 0) { + LOGP(fatal, "RICH quadrant geometry requires nTiles to be positive and divisible by four; received {}", richPars.nTiles); + } + if (vesselPhiGap < 0.0 || vesselThicknessShieldingLateral <= 0.0) { + LOGP(fatal, "RICH quadrant gap must be non-negative and lateral shielding thickness must be positive"); + } + const double totalBoundaryWidth = 2.0 * vesselThicknessShieldingLateral + vesselPhiGap; + if (totalBoundaryWidth >= 2.0 * richPars.rMin || vesselPhiGap >= 2.0 * shieldRMin) { + LOGP(fatal, "RICH quadrant boundary dimensions are incompatible with rMin={} cm and shieldRMin={} cm", richPars.rMin, shieldRMin); + } + const double quadrantDeltaPhi = solveQuadrantDeltaPhi(richPars.nTiles, richPars.rMin, totalBoundaryWidth); + if (!(quadrantDeltaPhi > 0.0)) { + LOGP(fatal, "RICH could not solve the quadrant module angular pitch"); + } + } + + // Name of the gas mother volume. This name will also be passed + // to each Ring so that the ring components become its daughters. + const char* richGasMotherName = "RICH_GAS_MOTHER"; + + TGeoMedium* medCO2 = gGeoManager->GetMedium("RCH_CO2$"); + if (!medCO2) { + LOGP(fatal, "RICH: CO2 medium not found"); + } + + TGeoMedium* medPeek = gGeoManager->GetMedium("RCH_PEEK$"); + if (!medPeek) { + LOGP(fatal, "RICH: PEEK medium not found"); + } + + TGeoMedium* medArmaFlex = gGeoManager->GetMedium("RCH_ARMAFLEX$"); + if (!medArmaFlex) { + LOGP(fatal, "RICH: ArmaFlex medium not found"); + } + + TGeoMedium* medArmaGel = gGeoManager->GetMedium("RCH_ARMAGEL$"); + if (!medArmaGel) { + LOGP(fatal, "RICH: ArmaGel medium not found"); + } + prepareLayout(); // Preparing the positions of the rings and tiles + // The gas mother includes the side-wall region and both end caps. ( as vessel ) + const double gasEnvelopeLengthZ = shieldLengthZ + 2.0 * endCapThicknessZ; + auto* gasEnvelopeShape = new TGeoTube(shieldRMin, shieldRMax, gasEnvelopeLengthZ / 2.0); + auto* gasEnvelopeVolume = new TGeoVolume(richGasMotherName, gasEnvelopeShape, medCO2); + + gasEnvelopeVolume->SetLineColor(kBlue - 9); + gasEnvelopeVolume->SetTransparency(90); + + // The gas envelope is a daughter of the general RICH volume. + vRICH->AddNode(gasEnvelopeVolume, 1, new TGeoTranslation(0.0, 0.0, 0.0)); + + if (!flagUseQuadrants) { + // ============================================================ + // Inner cylindrical insulating wall + // + // Radial interval: + // shieldRMin --> shieldRMin + innerWallThickness + // + // Longitudinal interval: + // -shieldLengthZ/2 --> +shieldLengthZ/2 + // ============================================================ + auto* innerWallShape = new TGeoTube(shieldRMin, shieldRMin + innerWallThickness, shieldLengthZ / 2.0); + auto* innerWallVolume = new TGeoVolume("RICH_SHIELD_INNER_WALL", innerWallShape, medArmaGel); + + innerWallVolume->SetLineColor(kOrange - 8); // kGray + innerWallVolume->SetTransparency(0); // 80 + gasEnvelopeVolume->AddNode(innerWallVolume, 1, new TGeoTranslation(0.0, 0.0, 0.0)); + + // ============================================================ + // Outer cylindrical insulating wall + // + // Radial interval: + // shieldRMax - outerWallThickness --> shieldRMax + // + // Longitudinal interval: + // -shieldLengthZ/2 --> +shieldLengthZ/2 + // ============================================================ + + auto* outerWallShape = new TGeoTube(shieldRMax - outerWallThickness, shieldRMax, shieldLengthZ / 2.0); + auto* outerWallVolume = new TGeoVolume("RICH_SHIELD_OUTER_WALL", outerWallShape, medArmaGel); + + outerWallVolume->SetLineColor(kOrange - 8); // kGray + outerWallVolume->SetTransparency(0); // 80 + gasEnvelopeVolume->AddNode(outerWallVolume, 1, new TGeoTranslation(0.0, 0.0, 0.0)); + + // ============================================================ + // Insulating end caps + // + // Each end cap covers: + // shieldRMin --> shieldRMax + // + // Each has full thickness: + // endCapThicknessZ + // ============================================================ + + auto* endCapShape = new TGeoTube(shieldRMin, shieldRMax, endCapThicknessZ / 2.0); + auto* endCapPlusVolume = new TGeoVolume("RICH_SHIELD_ENDCAP_PLUS", endCapShape, medArmaGel); + auto* endCapMinusVolume = new TGeoVolume("RICH_SHIELD_ENDCAP_MINUS", endCapShape, medArmaGel); + + endCapPlusVolume->SetLineColor(kOrange - 8); // kGray + endCapPlusVolume->SetTransparency(0); // 80 + + endCapMinusVolume->SetLineColor(kOrange - 8); // kGray + endCapMinusVolume->SetTransparency(0); // 80 + + const double endCapCenterZ = shieldLengthZ / 2.0 + endCapThicknessZ / 2.0; + + gasEnvelopeVolume->AddNode(endCapPlusVolume, 1, new TGeoTranslation(0.0, 0.0, endCapCenterZ)); + gasEnvelopeVolume->AddNode(endCapMinusVolume, 1, new TGeoTranslation(0.0, 0.0, -endCapCenterZ)); + } else { + // ============================================================ + // Four independent insulating vessel quadrants + // ============================================================ + const double totalBoundaryWidth = 2.0 * vesselThicknessShieldingLateral + vesselPhiGap; + const double moduleDeltaPhi = solveQuadrantDeltaPhi(richPars.nTiles, richPars.rMin, totalBoundaryWidth); + const double moduleExtraPhi = TMath::ASin(totalBoundaryWidth / (2.0 * richPars.rMin)); + const double vesselGapHalfPhi = TMath::ASin(vesselPhiGap / (2.0 * shieldRMin)); + const double quadrantSpanPhi = TMath::Pi() / 2.0 - 2.0 * vesselGapHalfPhi; + const int modulesPerQuadrant = richPars.nTiles / 4; + + // Remaining angular space between the last module of a quadrant and the following vessel gap. + const double endModuleExtraPhi = TMath::Pi() / 2.0 - moduleExtraPhi - static_cast(modulesPerQuadrant) * moduleDeltaPhi; + const double lateralStartWallSpanPhi = moduleExtraPhi - vesselGapHalfPhi; + const double lateralEndWallSpanPhi = endModuleExtraPhi - vesselGapHalfPhi; + + if (quadrantSpanPhi <= 0.0 || lateralStartWallSpanPhi <= 0.0 || lateralEndWallSpanPhi <= 0.0 || lateralStartWallSpanPhi + lateralEndWallSpanPhi >= quadrantSpanPhi) { + LOGP(fatal, "RICH invalid quadrant angular dimensions: vessel span={}, start wall span={}, end wall span={}", quadrantSpanPhi, lateralStartWallSpanPhi, lateralEndWallSpanPhi); + } + + const double radToDeg = 180.0 / TMath::Pi(); + const double quadrantSpanDeg = quadrantSpanPhi * radToDeg; + const double lateralStartWallSpanDeg = lateralStartWallSpanPhi * radToDeg; + const double lateralEndWallSpanDeg = lateralEndWallSpanPhi * radToDeg; + const double vesselGapHalfDeg = vesselGapHalfPhi * radToDeg; + const double innerGasRadius = shieldRMin + innerWallThickness; + const double outerGasRadius = shieldRMax - outerWallThickness; + + // Inner cylindrical shielding, divided into four sectors. + auto* innerWallQuadrantShape = new TGeoTubeSeg("RICH_SHIELD_INNER_WALL_QUADRANT_SHAPE", shieldRMin, innerGasRadius, shieldLengthZ / 2.0, 0.0, quadrantSpanDeg); + + // Outer cylindrical shielding, divided into four sectors. + auto* outerWallQuadrantShape = new TGeoTubeSeg("RICH_SHIELD_OUTER_WALL_QUADRANT_SHAPE", outerGasRadius, shieldRMax, shieldLengthZ / 2.0, 0.0, quadrantSpanDeg); + + // End caps divided into four sectors. + auto* endCapQuadrantShape = new TGeoTubeSeg("RICH_SHIELD_ENDCAP_QUADRANT_SHAPE", shieldRMin, shieldRMax, endCapThicknessZ / 2.0, 0.0, quadrantSpanDeg); + + // Lateral wall at the beginning of each quadrant. + auto* lateralStartWallShape = new TGeoTubeSeg("RICH_SHIELD_LATERAL_START_WALL_SHAPE", innerGasRadius, outerGasRadius, shieldLengthZ / 2.0, 0.0, lateralStartWallSpanDeg); + + // Lateral wall at the end of each quadrant. + auto* lateralEndWallShape = new TGeoTubeSeg("RICH_SHIELD_LATERAL_END_WALL_SHAPE", innerGasRadius, outerGasRadius, shieldLengthZ / 2.0, 0.0, lateralEndWallSpanDeg); + + auto* innerWallQuadrantVolume = new TGeoVolume("RICH_SHIELD_INNER_WALL_QUADRANT", innerWallQuadrantShape, medArmaGel); + auto* outerWallQuadrantVolume = new TGeoVolume("RICH_SHIELD_OUTER_WALL_QUADRANT", outerWallQuadrantShape, medArmaGel); + auto* endCapPlusQuadrantVolume = new TGeoVolume("RICH_SHIELD_ENDCAP_PLUS_QUADRANT", endCapQuadrantShape, medArmaGel); + auto* endCapMinusQuadrantVolume = new TGeoVolume("RICH_SHIELD_ENDCAP_MINUS_QUADRANT", endCapQuadrantShape, medArmaGel); + auto* lateralStartWallVolume = new TGeoVolume("RICH_SHIELD_LATERAL_START_WALL", lateralStartWallShape, medArmaGel); + auto* lateralEndWallVolume = new TGeoVolume("RICH_SHIELD_LATERAL_END_WALL", lateralEndWallShape, medArmaGel); + + innerWallQuadrantVolume->SetLineColor(kOrange - 8); // kGray + outerWallQuadrantVolume->SetLineColor(kOrange - 8); // kGray + endCapPlusQuadrantVolume->SetLineColor(kOrange - 8); // kGray + endCapMinusQuadrantVolume->SetLineColor(kOrange - 8); // kGray + lateralStartWallVolume->SetLineColor(kOrange - 8); // kGray + lateralEndWallVolume->SetLineColor(kOrange - 8); // kGray + + innerWallQuadrantVolume->SetTransparency(0); // 80 + outerWallQuadrantVolume->SetTransparency(0); // 80 + endCapPlusQuadrantVolume->SetTransparency(0); // 80 + endCapMinusQuadrantVolume->SetTransparency(0); // 80 + lateralStartWallVolume->SetTransparency(0); // 80 + lateralEndWallVolume->SetTransparency(0); // 80 + + const double endCapCenterZ = shieldLengthZ / 2.0 + endCapThicknessZ / 2.0; + + for (int quadrant = 0; quadrant < 4; quadrant++) { + + const double quadrantStartDeg = -45.0 + static_cast(quadrant) * 90.0 + vesselGapHalfDeg; + const double quadrantEndDeg = quadrantStartDeg + quadrantSpanDeg; + + auto makeRotation = [&](const char* prefix, double angleDeg) { + auto* rotation = new TGeoRotation(Form("%s_%d", prefix, quadrant)); + rotation->RotateZ(angleDeg); + return rotation; + }; + + // Inner cylindrical wall sector. + gasEnvelopeVolume->AddNode(innerWallQuadrantVolume, quadrant + 1, new TGeoCombiTrans(0.0, 0.0, 0.0, makeRotation("RICHInnerQuadrantRotation", quadrantStartDeg))); + // Outer cylindrical wall sector. + gasEnvelopeVolume->AddNode(outerWallQuadrantVolume, quadrant + 1, new TGeoCombiTrans(0.0, 0.0, 0.0, makeRotation("RICHOuterQuadrantRotation", quadrantStartDeg))); + // Positive-z end cap sector. + gasEnvelopeVolume->AddNode(endCapPlusQuadrantVolume, quadrant + 1, new TGeoCombiTrans(0.0, 0.0, endCapCenterZ, makeRotation("RICHEndCapPlusQuadrantRotation", quadrantStartDeg))); + // Negative-z end cap sector. + gasEnvelopeVolume->AddNode(endCapMinusQuadrantVolume, quadrant + 1, new TGeoCombiTrans(0.0, 0.0, -endCapCenterZ, makeRotation("RICHEndCapMinusQuadrantRotation", quadrantStartDeg))); + // Start-side lateral wall. + gasEnvelopeVolume->AddNode(lateralStartWallVolume, quadrant + 1, new TGeoCombiTrans(0.0, 0.0, 0.0, makeRotation("RICHLateralStartRotation", quadrantStartDeg))); + // End-side lateral wall. + gasEnvelopeVolume->AddNode(lateralEndWallVolume, quadrant + 1, new TGeoCombiTrans(0.0, 0.0, 0.0, makeRotation("RICHLateralEndRotation", quadrantEndDeg - lateralEndWallSpanDeg))); + } + + LOGP(info, "RICH quadrant geometry: module pitch={} deg, module boundary half-gap={} deg, vessel half-gap={} deg", moduleDeltaPhi * radToDeg, moduleExtraPhi * radToDeg, vesselGapHalfDeg); + } + + // ============================================================ modules for (int iRing{0}; iRing < richPars.nRings; ++iRing) { if (!richPars.oddGeom && iRing == (richPars.nRings / 2)) { continue; @@ -156,18 +669,18 @@ void Detector::createGeometry() richPars.rMin, richPars.rMax, richPars.radiatorThickness, - (float)mVTile1[iRing], - (float)mVTile2[iRing], - (float)mLAerogelZ[iRing], + (double)mVTile1[iRing], + (double)mVTile2[iRing], + (double)mLAerogelZ[iRing], richPars.detectorThickness, - (float)mVMirror1[iRing], - (float)mVMirror2[iRing], + (double)mVMirror1[iRing], + (double)mVMirror2[iRing], richPars.zBaseSize, - (float)mR0Radiator[iRing], - (float)mR0PhotoDet[iRing], - (float)mTRplusG[iRing], - (float)mThetaBi[iRing], - GeometryTGeo::getRICHVolPattern()}; + (double)mR0Radiator[iRing], + (double)mR0PhotoDet[iRing], + (double)mTRplusG[iRing], + (double)mThetaBi[iRing], + richGasMotherName}; // GeometryTGeo::getRICHVolPattern() } if (richPars.enableFWDRich) { @@ -233,7 +746,12 @@ void Detector::Reset() bool Detector::ProcessHits(FairVolume* vol) { // This method is called from the MC stepping - if (!(fMC->TrackCharge())) { + + constexpr int kOpticalPhotonPDG = 50000050; + const bool isOpticalPhoton = (fMC->TrackPid() == kOpticalPhotonPDG); + const bool isChargedParticle = (TMath::Abs(fMC->TrackCharge()) > 0.0); + // Reject neutral particles other than optical photons. + if (!isChargedParticle && !isOpticalPhoton) { return false; } @@ -242,6 +760,39 @@ bool Detector::ProcessHits(FairVolume* vol) // Is it needed to keep a track reference when the outer ITS volume is encountered? auto stack = (o2::data::Stack*)fMC->GetStack(); + + // Only the active silicon volumes are registered as sensitive in + // defineSensitiveVolumes(). The explicit volume-name check is kept + // as a safety guard in case additional sensitive volumes are added. + if (isOpticalPhoton) { + const char* currentVolumeName = fMC->CurrentVolName(); + const bool isActiveSiliconVolume = currentVolumeName && TString(currentVolumeName).BeginsWith(GeometryTGeo::getRICHSensorPattern()); + // Create only one hit when entering the active silicon. + if (!isActiveSiliconVolume || !fMC->IsTrackEntering()) { + return false; + } + TLorentzVector photonPosition; + TLorentzVector photonMomentum; + fMC->TrackPosition(photonPosition); + fMC->TrackMomentum(photonMomentum); + constexpr unsigned char photonStatus = Hit::kTrackEntering; + addHit( + stack->GetCurrentTrackNumber(), + lay, + photonPosition.Vect(), + photonPosition.Vect(), + photonMomentum.Vect(), + photonMomentum.E(), + photonPosition.T(), + 0.0, + photonStatus, + photonStatus); + + stack->addHit(GetDetId()); + + return true; + } + if (fMC->IsTrackExiting() && (lay == 0 || lay == mRings.size() - 1)) { // Keep the track refs for the innermost and outermost rings only o2::TrackReference tr(*fMC, GetDetId()); @@ -380,30 +931,121 @@ void Detector::prepareLayout() } // Dimensioning tiles - double percentage = 0.999; - for (int iRing = 0; iRing < richPars.nRings; iRing++) { - if (iRing == richPars.nRings / 2) { - mVMirror1[iRing] = percentage * 2.0 * richPars.rMax * TMath::Sin(TMath::Pi() / double(richPars.nTiles)); - mVMirror2[iRing] = percentage * 2.0 * richPars.rMax * TMath::Sin(TMath::Pi() / double(richPars.nTiles)); - mVTile1[iRing] = percentage * 2.0 * richPars.rMin * TMath::Sin(TMath::Pi() / double(richPars.nTiles)); - mVTile2[iRing] = percentage * 2.0 * richPars.rMin * TMath::Sin(TMath::Pi() / double(richPars.nTiles)); - } else if (iRing > richPars.nRings / 2) { - mVMirror1[iRing] = percentage * 2.0 * richPars.rMax * TMath::Sin(TMath::Pi() / double(richPars.nTiles)); - mVMirror2[iRing] = percentage * 2.0 * mMinRadialMirror[iRing] * TMath::Sin(TMath::Pi() / double(richPars.nTiles)); - mVTile1[iRing] = percentage * 2.0 * mMaxRadialRadiator[iRing] * TMath::Sin(TMath::Pi() / double(richPars.nTiles)); - mVTile2[iRing] = percentage * 2.0 * richPars.rMin * TMath::Sin(TMath::Pi() / double(richPars.nTiles)); - } else if (iRing < richPars.nRings / 2) { - mVMirror2[iRing] = percentage * 2.0 * richPars.rMax * TMath::Sin(TMath::Pi() / double(richPars.nTiles)); - mVMirror1[iRing] = percentage * 2.0 * mMinRadialMirror[iRing] * TMath::Sin(TMath::Pi() / double(richPars.nTiles)); - mVTile2[iRing] = percentage * 2.0 * mMaxRadialRadiator[iRing] * TMath::Sin(TMath::Pi() / double(richPars.nTiles)); - mVTile1[iRing] = percentage * 2.0 * richPars.rMin * TMath::Sin(TMath::Pi() / double(richPars.nTiles)); + if (!richPars.flagUseQuadrants) { + double percentage = 0.999; + for (int iRing = 0; iRing < richPars.nRings; iRing++) { + if (iRing == richPars.nRings / 2) { + mVMirror1[iRing] = percentage * 2.0 * richPars.rMax * TMath::Tan(TMath::Pi() / double(richPars.nTiles)); + mVMirror2[iRing] = percentage * 2.0 * richPars.rMax * TMath::Tan(TMath::Pi() / double(richPars.nTiles)); + mVTile1[iRing] = percentage * 2.0 * richPars.rMin * TMath::Tan(TMath::Pi() / double(richPars.nTiles)); + mVTile2[iRing] = percentage * 2.0 * richPars.rMin * TMath::Tan(TMath::Pi() / double(richPars.nTiles)); + } else if (iRing > richPars.nRings / 2) { + mVMirror1[iRing] = percentage * 2.0 * richPars.rMax * TMath::Tan(TMath::Pi() / double(richPars.nTiles)); + mVMirror2[iRing] = percentage * 2.0 * mMinRadialMirror[iRing] * TMath::Tan(TMath::Pi() / double(richPars.nTiles)); + mVTile1[iRing] = percentage * 2.0 * mMaxRadialRadiator[iRing] * TMath::Tan(TMath::Pi() / double(richPars.nTiles)); + mVTile2[iRing] = percentage * 2.0 * richPars.rMin * TMath::Tan(TMath::Pi() / double(richPars.nTiles)); + } else { + mVMirror2[iRing] = percentage * 2.0 * richPars.rMax * TMath::Tan(TMath::Pi() / double(richPars.nTiles)); + mVMirror1[iRing] = percentage * 2.0 * mMinRadialMirror[iRing] * TMath::Tan(TMath::Pi() / double(richPars.nTiles)); + mVTile2[iRing] = percentage * 2.0 * mMaxRadialRadiator[iRing] * TMath::Tan(TMath::Pi() / double(richPars.nTiles)); + mVTile1[iRing] = percentage * 2.0 * richPars.rMin * TMath::Tan(TMath::Pi() / double(richPars.nTiles)); + } + } + + } else { + + const double totalBoundaryWidth = 2.0 * richPars.vesselThicknessShieldingLateral + richPars.vesselPhiGap; + const double quadrantDeltaPhi = solveQuadrantDeltaPhi(richPars.nTiles, richPars.rMin, totalBoundaryWidth); + if (!(quadrantDeltaPhi > 0.0)) { + LOGP(fatal, "RICH could not solve the quadrant module angular pitch"); + } + const double halfWidthFactor = TMath::Tan(quadrantDeltaPhi / 2.0); + double percentage = 0.999; + for (int iRing = 0; iRing < richPars.nRings; iRing++) { + if (iRing == richPars.nRings / 2) { + mVMirror1[iRing] = percentage * 2.0 * richPars.rMax * halfWidthFactor; + mVMirror2[iRing] = percentage * 2.0 * richPars.rMax * halfWidthFactor; + mVTile1[iRing] = percentage * 2.0 * richPars.rMin * halfWidthFactor; + mVTile2[iRing] = percentage * 2.0 * richPars.rMin * halfWidthFactor; + } else if (iRing > richPars.nRings / 2) { + mVMirror1[iRing] = percentage * 2.0 * richPars.rMax * halfWidthFactor; + mVMirror2[iRing] = percentage * 2.0 * mMinRadialMirror[iRing] * halfWidthFactor; + mVTile1[iRing] = percentage * 2.0 * mMaxRadialRadiator[iRing] * halfWidthFactor; + mVTile2[iRing] = percentage * 2.0 * richPars.rMin * halfWidthFactor; + + } else { + mVMirror2[iRing] = percentage * 2.0 * richPars.rMax * halfWidthFactor; + mVMirror1[iRing] = percentage * 2.0 * mMinRadialMirror[iRing] * halfWidthFactor; + mVTile2[iRing] = percentage * 2.0 * mMaxRadialRadiator[iRing] * halfWidthFactor; + mVTile1[iRing] = percentage * 2.0 * richPars.rMin * halfWidthFactor; + } + } + } + + // ============================================================ + // Cylindrical aerogel geometry + // ============================================================ + // + // In this mode the photosensors remain projective, but all + // aerogel tiles: + // + // - have identical dimensions; + // - are parallel to the beam axis; + // - lie at the same cylindrical radius; + // - are uniformly distributed along Z. + // + if (richPars.useCylindricalAerogel) { + + // In the even geometry the central projective ring is skipped + // in createGeometry(), so the number of actual aerogel rows is + // nRings - 1. + const int nAerogelRows = richPars.oddGeom ? richPars.nRings : richPars.nRings - 1; + + if (nAerogelRows <= 0) { + LOGP(fatal, "Invalid number of cylindrical aerogel rows: {}", nAerogelRows); + } + + if (richPars.nTiles <= 0) { + LOGP(fatal, "Invalid number of aerogel tiles in phi: {}", richPars.nTiles); + } + + const double thetaRef = 2.0 * TMath::ATan(TMath::Exp(-richPars.cylindricalAerogelEtaRef)); + const double cylindricalAerogelTileSizeZ = (2.0 * richPars.rMin / TMath::Tan(thetaRef)) / static_cast(nAerogelRows); + // const double cylindricalAerogelTileSizeRPhi = 2.0 * richPars.rMin * TMath::Tan(TMath::Pi() / static_cast(richPars.nTiles)); + double cylindricalAerogelTileSizeRPhi = 0.0; + + if (!richPars.flagUseQuadrants) { + // Original uniform-phi geometry. + cylindricalAerogelTileSizeRPhi = 2.0 * richPars.rMin * TMath::Tan(TMath::Pi() / static_cast(richPars.nTiles)); + } else { + const double totalBoundaryWidth = 2.0 * richPars.vesselThicknessShieldingLateral + richPars.vesselPhiGap; + const double quadrantDeltaPhi = solveQuadrantDeltaPhi(richPars.nTiles, richPars.rMin, totalBoundaryWidth); + cylindricalAerogelTileSizeRPhi = 2.0 * richPars.rMin * TMath::Tan(quadrantDeltaPhi / 2.0); + } + + LOGP(info, "Cylindrical aerogel: rows={}, etaRef={}, tileSizeZ={} cm, tileSizeRPhi={} cm", nAerogelRows, richPars.cylindricalAerogelEtaRef, cylindricalAerogelTileSizeZ, cylindricalAerogelTileSizeRPhi); + + for (int iRing = 0; iRing < richPars.nRings; iRing++) { + mLAerogelZ[iRing] = cylindricalAerogelTileSizeZ; + + // Equal values make the TGeoArb8 a rectangle instead of the projective trapezoid. + mVTile1[iRing] = cylindricalAerogelTileSizeRPhi; + mVTile2[iRing] = cylindricalAerogelTileSizeRPhi; } } // Translation parameters for (size_t iRing{0}; iRing < richPars.nRings; ++iRing) { - mR0Radiator[iRing] = mR0Tilt[iRing] - (mTRplusG[iRing] - richPars.radiatorThickness / 2) * TMath::Cos(mThetaBi[iRing]); - mR0PhotoDet[iRing] = mR0Tilt[iRing] - (richPars.detectorThickness / 2) * TMath::Cos(mThetaBi[iRing]); + + if (richPars.useCylindricalAerogel) { + mR0Radiator[iRing] = richPars.rMin + richPars.radiatorThickness / 2.0; + } else { + // Original projective aerogel position. + mR0Radiator[iRing] = mR0Tilt[iRing] - (mTRplusG[iRing] - richPars.radiatorThickness / 2.0) * TMath::Cos(mThetaBi[iRing]); + } + + // Photosensors remain projective for both configurations. + mR0PhotoDet[iRing] = mR0Tilt[iRing] - richPars.detectorThickness / 2.0 * TMath::Cos(mThetaBi[iRing]); } // FWD and BWD RICH diff --git a/Detectors/Upgrades/ALICE3/RICH/simulation/src/RICHRing.cxx b/Detectors/Upgrades/ALICE3/RICH/simulation/src/RICHRing.cxx index 46e6ae515b390..45c5c52b27043 100644 --- a/Detectors/Upgrades/ALICE3/RICH/simulation/src/RICHRing.cxx +++ b/Detectors/Upgrades/ALICE3/RICH/simulation/src/RICHRing.cxx @@ -19,183 +19,594 @@ #include #include +#include +#include + namespace o2 { namespace rich { +namespace // quadrant operations +{ + +double quadrantDeltaPhiEquation(double x, int nTilesPhi, double rMin, double totalBoundaryWidth) +{ + const double argument = totalBoundaryWidth * TMath::Cos(x / 2.0) / (2.0 * rMin); + if (TMath::Abs(argument) >= 1.0) { + return std::numeric_limits::quiet_NaN(); + } + const double rhs = 2.0 * TMath::Pi() / static_cast(nTilesPhi) - (8.0 / static_cast(nTilesPhi)) * TMath::ASin(argument); + return rhs - x; +} + +double solveQuadrantDeltaPhi(int nTilesPhi, double rMin, double totalBoundaryWidth) +{ + double lower = 0.0; + double upper = 1.1 * 2.0 * TMath::Pi() / static_cast(nTilesPhi); + double fLower = quadrantDeltaPhiEquation(lower, nTilesPhi, rMin, totalBoundaryWidth); + double fUpper = quadrantDeltaPhiEquation(upper, nTilesPhi, rMin, totalBoundaryWidth); + if (!std::isfinite(fLower) || !std::isfinite(fUpper) || fLower * fUpper > 0.0) { + return -1.0; + } + constexpr double tolerance = 1.0e-12; + constexpr int maxIterations = 200; + for (int iteration = 0; iteration < maxIterations; iteration++) { + const double middle = 0.5 * (lower + upper); + const double fMiddle = quadrantDeltaPhiEquation(middle, nTilesPhi, rMin, totalBoundaryWidth); + if (!std::isfinite(fMiddle)) { + return -1.0; + } + if (TMath::Abs(fMiddle) < tolerance || 0.5 * (upper - lower) < tolerance) { + return middle; + } + if (fLower * fMiddle < 0.0) { + upper = middle; + fUpper = fMiddle; + } else { + lower = middle; + fLower = fMiddle; + } + } + return 0.5 * (lower + upper); +} + +double quadrantModulePhi(int moduleIndex, int nTilesPhi, double deltaPhi, double extraPhi) +{ + const int modulesPerQuadrant = nTilesPhi / 4; + const int quadrant = moduleIndex / modulesPerQuadrant; + return extraPhi + static_cast(moduleIndex) * deltaPhi - TMath::Pi() / 4.0 + deltaPhi / 2.0 + 2.0 * static_cast(quadrant) * extraPhi; +} + +} // namespace + Ring::Ring(int rPosId, int nTilesPhi, - float rMin, - float rMax, - float radThick, - float radYmin, - float radYmax, - float radZ, - float photThick, - float photYmin, - float photYmax, - float photZ, - float radRad0, - float photR0, - float aerDetDistance, - float thetaB, + double rMin, + double rMax, + double radThick, + double radYmin, + double radYmax, + double radZ, + double photThick, + double photYmin, + double photYmax, + double photZ, + double radRad0, + double photR0, + double aerDetDistance, + double thetaB, const std::string motherName) : mNTiles{nTilesPhi}, mPosId{rPosId}, mRadThickness{radThick} { TGeoManager* geoManager = gGeoManager; TGeoVolume* motherVolume = geoManager->GetVolume(motherName.c_str()); + + if (!motherVolume) { + LOGP(fatal, + "RICH: mother volume {} not found while creating ring {}", + motherName, + rPosId); + } + + const auto& richPars = RICHBaseParam::Instance(); + + const bool useCylindricalAerogel = richPars.useCylindricalAerogel; + TGeoMedium* medAerogel = gGeoManager->GetMedium("RCH_AEROGEL$"); if (!medAerogel) { LOGP(fatal, "RICH: Aerogel medium not found"); } + TGeoMedium* medSi = gGeoManager->GetMedium("RCH_SILICON$"); if (!medSi) { LOGP(fatal, "RICH: Silicon medium not found"); } + + TGeoMedium* medCO2 = gGeoManager->GetMedium("RCH_CO2$"); + if (!medCO2) { + LOGP(fatal, "RICH: CO2 medium not found"); + } + + TGeoMedium* medFR4 = gGeoManager->GetMedium("RCH_FR4$"); + if (!medFR4) { + LOGP(fatal, "RICH: FR4 medium not found"); + } + TGeoMedium* medAr = gGeoManager->GetMedium("RCH_ARGON$"); if (!medAr) { LOGP(fatal, "RICH: Argon medium not found"); } - std::vector radiatorTiles(nTilesPhi), photoTiles(nTilesPhi), argonSectors(nTilesPhi); + + TGeoMedium* medAl = gGeoManager->GetMedium("RCH_ALUMINUM$"); + if (!medAl) { + LOGP(fatal, "RICH: Aluminum medium not found"); + } + + TGeoMedium* medSiAbsorber = gGeoManager->GetMedium("RCH_SILICON_ABSORBER$"); + if (!medSiAbsorber) { + LOGP(fatal, "RICH: Passive silicon absorber medium not found"); + } + + TGeoMedium* medSilicone = gGeoManager->GetMedium("RCH_SILICONE$"); + if (!medSilicone) { + LOGP(fatal, "RICH: Silicone medium not found"); + } + + TGeoMedium* medHTCC = gGeoManager->GetMedium("RCH_HTCC$"); + if (!medHTCC) { + LOGP(fatal, "RICH: HTCC medium not found"); + } + + std::vector radiatorTiles(nTilesPhi), photoFrames(nTilesPhi), photoTiles(nTilesPhi), gasSectors(nTilesPhi); LOGP(info, "Creating ring: id: {} with {} tiles. ", rPosId, nTilesPhi); LOGP(info, "Rmin: {} Rmax: {} RadThick: {} RadYmin: {} RadYmax: {} RadZ: {} PhotThick: {} PhotYmin: {} PhotYmax: {} PhotZ: {}, zTransRad: {}, zTransPhot: {}, ThetaB: {}", rMin, rMax, radThick, radYmin, radYmax, radZ, photThick, photYmin, photYmax, photZ, radRad0, photR0, thetaB); - float deltaPhiDeg = 360.0 / nTilesPhi; // Transformation are constructed in degrees... - float thetaBDeg = thetaB * 180.0 / TMath::Pi(); - int radTileCount{0}, photTileCount{0}, argSectorsCount{0}; + // Use different phi depending on use of quadrants or not + const bool flagUseQuadrants = richPars.flagUseQuadrants; + if (flagUseQuadrants && (nTilesPhi <= 0 || nTilesPhi % 4 != 0)) { + LOGP(fatal, "RICH quadrant geometry requires nTilesPhi to be positive and divisible by four; received {}", nTilesPhi); + } + const double regularDeltaPhi = 2.0 * TMath::Pi() / static_cast(nTilesPhi); + double moduleDeltaPhi = regularDeltaPhi; + double quadrantExtraPhi = 0.0; + + if (flagUseQuadrants) { + const double totalBoundaryWidth = 2.0 * richPars.vesselThicknessShieldingLateral + richPars.vesselPhiGap; + if (totalBoundaryWidth >= 2.0 * rMin) { + LOGP(fatal, "RICH quadrant boundary width {} cm is incompatible with rMin={} cm", totalBoundaryWidth, rMin); + } + moduleDeltaPhi = solveQuadrantDeltaPhi(nTilesPhi, rMin, totalBoundaryWidth); + + quadrantExtraPhi = TMath::ASin(totalBoundaryWidth / (2.0 * rMin)); + + if (!(moduleDeltaPhi > 0.0)) { + LOGP(fatal, "RICH ring {} could not solve the quadrant angular pitch", rPosId); + } + } + + auto modulePhiRad = [&](int moduleIndex) { + if (!flagUseQuadrants) { + // Original placement exactly. + return static_cast(moduleIndex) * regularDeltaPhi; + } + return quadrantModulePhi(moduleIndex, nTilesPhi, moduleDeltaPhi, quadrantExtraPhi); + }; + + const double thetaBDeg = thetaB * 180.0 / TMath::Pi(); + + const double sipmActiveSizeZ = richPars.sipmActiveSizeZ; + // const double sipmActiveSizeRPhi = richPars.sipmActiveSizeRPhi; + // Select width depending on having quadrants or not (and wall thickness) + const double sipmActiveSizeRPhi = flagUseQuadrants ? richPars.quadrantModuleSizeRPhi : richPars.sipmActiveSizeRPhi; + + const double pcb1Thickness = richPars.pcb1Thickness; + const double coolingPlateThickness = richPars.coolingPlateThickness; + const double pcb2Thickness = richPars.pcb2Thickness; + const double pcb3Thickness = richPars.pcb3Thickness; + + const double gapSiPMToPCB1 = richPars.gapSiPMToPCB1; + const double gapPCB1ToCoolingPlate = richPars.gapPCB1ToCoolingPlate; + const double gapCoolingPlateToPCB2 = richPars.gapCoolingPlateToPCB2; + const double gapPCB2ToPCB3 = richPars.gapPCB2ToPCB3; + + const bool oddGeom = richPars.oddGeom; + const bool useRectangularModules = richPars.useRectangularModules; + + const int nRings = richPars.nRings; + + const double moduleClearanceZ = richPars.moduleClearanceZ; + const double moduleClearanceRPhi = richPars.moduleClearanceRPhi; + + const double siliconeLayerThickness = richPars.siliconeLayerThickness; + const double activeSiliconThickness = richPars.activeSiliconThickness; + const double passiveSiliconThickness = photThick - activeSiliconThickness; + + const double siliconFrontSurfaceOffset = -photThick / 2.0; + const double siliconeCenterOffset = siliconFrontSurfaceOffset - siliconeLayerThickness / 2.0; + const double activeSiliconCenterOffset = siliconFrontSurfaceOffset + activeSiliconThickness / 2.0; + const double passiveSiliconCenterOffset = siliconFrontSurfaceOffset + activeSiliconThickness + passiveSiliconThickness / 2.0; + + if (siliconeLayerThickness <= 0.0) { + LOGP(fatal, "RICH: siliconeLayerThickness must be positive"); + } + + if (activeSiliconThickness <= 0.0 || activeSiliconThickness >= photThick) { + LOGP(fatal, "RICH: activeSiliconThickness={} cm must be larger than zero and smaller than detectorThickness={} cm", activeSiliconThickness, photThick); + } + + if (passiveSiliconThickness <= 0.0) { + LOGP(fatal, "RICH: passive silicon thickness must be positive"); + } + + if (moduleClearanceZ < 0.0 || moduleClearanceRPhi < 0.0) { + LOGP(fatal, "RICH: module clearances cannot be negative"); + } + + if (photThick <= 0.0 || sipmActiveSizeZ <= 0.0 || sipmActiveSizeRPhi <= 0.0 || pcb1Thickness <= 0.0 || coolingPlateThickness <= 0.0 || pcb2Thickness <= 0.0 || pcb3Thickness <= 0.0) { + LOGP(fatal, "RICH: SiPM and readout-stack dimensions must be positive"); + } + + if (gapSiPMToPCB1 < 0.0 || gapPCB1ToCoolingPlate < 0.0 || gapCoolingPlateToPCB2 < 0.0 || gapPCB2ToPCB3 < 0.0) { + LOGP(fatal, "RICH: readout-stack gaps cannot be negative"); + } + + const double minimumFrameSizeRPhi = photYmin < photYmax ? photYmin : photYmax; + + if (sipmActiveSizeZ > photZ || sipmActiveSizeRPhi > minimumFrameSizeRPhi) { + LOGP(fatal, + "RICH: rectangular module {} x {} cm2 does not fit inside the trapezoidal sector {} x [{}, {}] cm2 for ring {}. " + "For quadrant mode reduce: quadrantModuleSizeRPhi.", + sipmActiveSizeZ, sipmActiveSizeRPhi, photZ, photYmin, photYmax, rPosId); + } + + // Number of actual aerogel rows. + const int nAerogelRows = oddGeom ? nRings : nRings - 1; + + // Convert the projective-ring ID into a contiguous aerogel-row index. + // Example for nRings=11 and even geometry: + // projective IDs: 0 1 2 3 4 [5 skipped] 6 7 8 9 10 + // aerogel index: 0 1 2 3 4 5 6 7 8 9 + int aerogelRowIndex = rPosId; + if (!oddGeom && rPosId > nRings / 2) { + --aerogelRowIndex; + } + + const double cylindricalAerogelCenterZ = -0.5 * static_cast(nAerogelRows) * radZ + 0.5 * radZ + static_cast(aerogelRowIndex) * radZ; + + int radTileCount{0}, photTileCount{0}; // argSectorsCount{0}; + + if (flagUseQuadrants) { + LOGP(info, "RICH ring {} quadrant placement: deltaPhi={} deg, boundary half-gap={} deg", rPosId, moduleDeltaPhi * 180.0 / TMath::Pi(), quadrantExtraPhi * 180.0 / TMath::Pi()); + } + // Radiator tiles for (auto& radiatorTile : radiatorTiles) { // Local Z is the thin (radial) dimension, looking outward from the IP // (previously this was local X, while for running with ACTS we need local Z). // The placement rotation below is adjusted by +90 deg about Y // to keep the tile in the same physical position. - radiatorTile = new TGeoArb8(radThick / 2); - radiatorTile->SetVertex(0, radZ / 2, -radYmin / 2); - radiatorTile->SetVertex(1, -radZ / 2, -radYmax / 2); - radiatorTile->SetVertex(2, -radZ / 2, radYmax / 2); - radiatorTile->SetVertex(3, radZ / 2, radYmin / 2); - radiatorTile->SetVertex(4, radZ / 2, -radYmin / 2); - radiatorTile->SetVertex(5, -radZ / 2, -radYmax / 2); - radiatorTile->SetVertex(6, -radZ / 2, radYmax / 2); - radiatorTile->SetVertex(7, radZ / 2, radYmin / 2); + if (useCylindricalAerogel) { + // Including gab between adjacent aerogel tiles + const double cylindricalTileSizeZ = radZ - moduleClearanceZ; + const double cylindricalTileYmin = radYmin - moduleClearanceRPhi; + const double cylindricalTileYmax = radYmax - moduleClearanceRPhi; + if (cylindricalTileSizeZ <= 0.0 || cylindricalTileYmin <= 0.0 || cylindricalTileYmax <= 0.0) { + LOGP(fatal, "RICH: cylindrical-aerogel clearances are larger than the tile dimensions for ring {}", rPosId); + } + radiatorTile = new TGeoArb8(radThick / 2); + radiatorTile->SetVertex(0, cylindricalTileSizeZ / 2, -cylindricalTileYmin / 2); + radiatorTile->SetVertex(1, -cylindricalTileSizeZ / 2, -cylindricalTileYmax / 2); + radiatorTile->SetVertex(2, -cylindricalTileSizeZ / 2, cylindricalTileYmax / 2); + radiatorTile->SetVertex(3, cylindricalTileSizeZ / 2, cylindricalTileYmin / 2); + radiatorTile->SetVertex(4, cylindricalTileSizeZ / 2, -cylindricalTileYmin / 2); + radiatorTile->SetVertex(5, -cylindricalTileSizeZ / 2, -cylindricalTileYmax / 2); + radiatorTile->SetVertex(6, -cylindricalTileSizeZ / 2, cylindricalTileYmax / 2); + radiatorTile->SetVertex(7, cylindricalTileSizeZ / 2, cylindricalTileYmin / 2); + } else { + // Original non-cylindrical tile definition. + radiatorTile = new TGeoArb8(radThick / 2); + radiatorTile->SetVertex(0, radZ / 2, -radYmin / 2); + radiatorTile->SetVertex(1, -radZ / 2, -radYmax / 2); + radiatorTile->SetVertex(2, -radZ / 2, radYmax / 2); + radiatorTile->SetVertex(3, radZ / 2, radYmin / 2); + radiatorTile->SetVertex(4, radZ / 2, -radYmin / 2); + radiatorTile->SetVertex(5, -radZ / 2, -radYmax / 2); + radiatorTile->SetVertex(6, -radZ / 2, radYmax / 2); + radiatorTile->SetVertex(7, radZ / 2, radYmin / 2); + } TGeoVolume* radiatorTileVol = new TGeoVolume(Form("radTile_%d_%d", rPosId, radTileCount), radiatorTile, medAerogel); - radiatorTileVol->SetLineColor(kOrange - 8); + radiatorTileVol->SetLineColor(kBlue - 9); radiatorTileVol->SetLineWidth(1); + // const double phiDeg = static_cast(radTileCount) * deltaPhiDeg; + // const double phiRad = static_cast(radTileCount) * 2.0 * TMath::Pi() / static_cast(nTilesPhi); + + const double phiRad = modulePhiRad(radTileCount); + const double phiDeg = phiRad * 180.0 / TMath::Pi(); + auto* rotRadiator = new TGeoRotation(Form("radTileRotation_%d_%d", radTileCount, rPosId)); - rotRadiator->RotateY(90.0 - thetaBDeg); // +90 compensates the X->Z swap of the tile's local axes - rotRadiator->RotateZ(radTileCount * deltaPhiDeg); - auto* rotTransRadiator = new TGeoCombiTrans(radRad0 * TMath::Cos(radTileCount * TMath::Pi() / (nTilesPhi / 2)), - radRad0 * TMath::Sin(radTileCount * TMath::Pi() / (nTilesPhi / 2)), - radRad0 * TMath::Tan(thetaB), - rotRadiator); + if (useCylindricalAerogel) { + // The TGeoArb8 local Z axis is the thin direction. + // RotateY(90 degrees) maps that thin local Z direction onto the global radial direction at phi=0. + // There is no thetaB tilt because the cylindrical aerogel tiles are parallel to the beam axis. + rotRadiator->RotateY(90.0); + } else { + // Original projective rotation. + rotRadiator->RotateY(90.0 - thetaBDeg); + } + + // Rotate the radial tile around the beam axis to its phi sector. + rotRadiator->RotateZ(phiDeg); + + const double radiatorCenterZ = useCylindricalAerogel ? cylindricalAerogelCenterZ : radRad0 * TMath::Tan(thetaB); + + auto* rotTransRadiator = new TGeoCombiTrans(radRad0 * TMath::Cos(phiRad), radRad0 * TMath::Sin(phiRad), radiatorCenterZ, rotRadiator); motherVolume->AddNode(radiatorTileVol, 1, rotTransRadiator); radTileCount++; } - // Photosensor tiles - for (auto& photoTile : photoTiles) { - // Local Z is the thin (radial) dimension, looking outward from the IP - photoTile = new TGeoArb8(photThick / 2); - photoTile->SetVertex(0, photZ / 2, -photYmin / 2); - photoTile->SetVertex(1, -photZ / 2, -photYmax / 2); - photoTile->SetVertex(2, -photZ / 2, photYmax / 2); - photoTile->SetVertex(3, photZ / 2, photYmin / 2); - photoTile->SetVertex(4, photZ / 2, -photYmin / 2); - photoTile->SetVertex(5, -photZ / 2, -photYmax / 2); - photoTile->SetVertex(6, -photZ / 2, photYmax / 2); - photoTile->SetVertex(7, photZ / 2, photYmin / 2); - - TGeoVolume* photoTileVol = new TGeoVolume(Form("%s_%d_%d", GeometryTGeo::getRICHSensorPattern(), rPosId, photTileCount), photoTile, medSi); - photoTileVol->SetLineColor(kOrange - 8); - photoTileVol->SetLineWidth(1); - - auto* rotPhoto = new TGeoRotation(Form("photoTileRotation_%d_%d", photTileCount, rPosId)); - rotPhoto->RotateY(90.0 - thetaBDeg); // +90 compensates the X->Z swap of the tile's local axes - rotPhoto->RotateZ(photTileCount * deltaPhiDeg); - auto* rotTransPhoto = new TGeoCombiTrans(photR0 * TMath::Cos(photTileCount * TMath::Pi() / (nTilesPhi / 2)), - photR0 * TMath::Sin(photTileCount * TMath::Pi() / (nTilesPhi / 2)), - photR0 * TMath::Tan(thetaB), - rotPhoto); - - motherVolume->AddNode(photoTileVol, 1, rotTransPhoto); - photTileCount++; - } - - // Argon sectors "connect" radiator and photosensor tiles, they are not really physical - for (auto& argonSector : argonSectors) { - float separation{(aerDetDistance - radThick - photThick)}; + // Photosensor tiles: legacy trapezoidal modules and rectangular modules + if (!useRectangularModules) { + for (auto& photoTile : photoTiles) { + const double phiRad = modulePhiRad(photTileCount); + const double phiDeg = phiRad * 180.0 / TMath::Pi(); + // Local Z is the thin (radial) dimension, looking outward from the IP + photoTile = new TGeoArb8(photThick / 2); + photoTile->SetVertex(0, photZ / 2, -photYmin / 2); + photoTile->SetVertex(1, -photZ / 2, -photYmax / 2); + photoTile->SetVertex(2, -photZ / 2, photYmax / 2); + photoTile->SetVertex(3, photZ / 2, photYmin / 2); + photoTile->SetVertex(4, photZ / 2, -photYmin / 2); + photoTile->SetVertex(5, -photZ / 2, -photYmax / 2); + photoTile->SetVertex(6, -photZ / 2, photYmax / 2); + photoTile->SetVertex(7, photZ / 2, photYmin / 2); + + TGeoVolume* photoTileVol = new TGeoVolume(Form("%s_%d_%d", GeometryTGeo::getRICHSensorPattern(), rPosId, photTileCount), photoTile, medSi); + photoTileVol->SetLineColor(kOrange + 2); + photoTileVol->SetLineWidth(1); + + auto* rotPhoto = new TGeoRotation(Form("photoTileRotation_%d_%d", photTileCount, rPosId)); + rotPhoto->RotateY(90.0 - thetaBDeg); // +90 compensates the X->Z swap of the tile's local axes + // rotPhoto->RotateZ(photTileCount * deltaPhiDeg); + rotPhoto->RotateZ(phiDeg); + // auto* rotTransPhoto = new TGeoCombiTrans(photR0 * TMath::Cos(photTileCount * TMath::Pi() / (nTilesPhi / 2)), + // photR0 * TMath::Sin(photTileCount * TMath::Pi() / (nTilesPhi / 2)), + // photR0 * TMath::Tan(thetaB), + // rotPhoto); + auto* rotTransPhoto = new TGeoCombiTrans(photR0 * TMath::Cos(phiRad), photR0 * TMath::Sin(phiRad), photR0 * TMath::Tan(thetaB), rotPhoto); + + motherVolume->AddNode(photoTileVol, 1, rotTransPhoto); + photTileCount++; + } + } else // <-- New gemetry with rectangular modules + { + // Photosensor tiles and readout stack + for (auto& photoTile : photoTiles) { + // const double phiDeg = static_cast(photTileCount) * deltaPhiDeg; + // const double phiRad = static_cast(photTileCount) * 2.0 * TMath::Pi() / static_cast(nTilesPhi); + const double phiRad = modulePhiRad(photTileCount); + const double phiDeg = phiRad * 180.0 / TMath::Pi(); + + const double photoCenterR = photR0; + const double photoCenterZ = photR0 * TMath::Tan(thetaB); + + // Unit vector normal to the projective plane, pointing away from the IP. Positive offset places layer behind the SiPM. + const double normalRadial = TMath::Cos(thetaB); + const double normalZ = TMath::Sin(thetaB); + + auto makeProjectiveRotation = [&](const char* prefix) { + auto* rotation = new TGeoRotation(Form("%sRotation_%d_%d", prefix, photTileCount, rPosId)); + rotation->RotateY(90.0 - thetaBDeg); // same orientation as the original photosensor + rotation->RotateZ(phiDeg); + return rotation; + }; + + const double frameSizeZ = photZ - moduleClearanceZ; + const double frameYmin = photYmin - moduleClearanceRPhi; + const double frameYmax = photYmax - moduleClearanceRPhi; + // Footprint of the frames with the configured clearances for overlaps + auto makeFrameFootprint = [&](double thickness) { + auto* shape = new TGeoArb8(thickness / 2.0); + shape->SetVertex(0, frameSizeZ / 2.0, -frameYmin / 2.0); + shape->SetVertex(1, -frameSizeZ / 2.0, -frameYmax / 2.0); + shape->SetVertex(2, -frameSizeZ / 2.0, frameYmax / 2.0); + shape->SetVertex(3, frameSizeZ / 2.0, frameYmin / 2.0); + shape->SetVertex(4, frameSizeZ / 2.0, -frameYmin / 2.0); + shape->SetVertex(5, -frameSizeZ / 2.0, -frameYmax / 2.0); + shape->SetVertex(6, -frameSizeZ / 2.0, frameYmax / 2.0); + shape->SetVertex(7, frameSizeZ / 2.0, frameYmin / 2.0); + return shape; + }; + + auto makeRectangularFootprint = [&](double thickness) { + auto* shape = new TGeoArb8(thickness / 2.0); + shape->SetVertex(0, sipmActiveSizeZ / 2.0, -sipmActiveSizeRPhi / 2.0); + shape->SetVertex(1, -sipmActiveSizeZ / 2.0, -sipmActiveSizeRPhi / 2.0); + shape->SetVertex(2, -sipmActiveSizeZ / 2.0, sipmActiveSizeRPhi / 2.0); + shape->SetVertex(3, sipmActiveSizeZ / 2.0, sipmActiveSizeRPhi / 2.0); + shape->SetVertex(4, sipmActiveSizeZ / 2.0, -sipmActiveSizeRPhi / 2.0); + shape->SetVertex(5, -sipmActiveSizeZ / 2.0, -sipmActiveSizeRPhi / 2.0); + shape->SetVertex(6, -sipmActiveSizeZ / 2.0, sipmActiveSizeRPhi / 2.0); + shape->SetVertex(7, sipmActiveSizeZ / 2.0, sipmActiveSizeRPhi / 2.0); + return shape; + }; + + auto addReadoutLayer = [&](const char* prefix, + double thickness, + double centerOffset, + TGeoMedium* medium, + Color_t lineColor, + bool useRectangularFootprint) { + auto* shape = useRectangularFootprint ? makeRectangularFootprint(thickness) : makeFrameFootprint(thickness); + auto* volume = new TGeoVolume(Form("%s_%d_%d", prefix, rPosId, photTileCount), shape, medium); + volume->SetLineColor(lineColor); + volume->SetLineWidth(1); + const double layerCenterR = photoCenterR + centerOffset * normalRadial; + const double layerCenterZ = photoCenterZ + centerOffset * normalZ; + auto* transform = new TGeoCombiTrans(layerCenterR * TMath::Cos(phiRad), layerCenterR * TMath::Sin(phiRad), layerCenterZ, makeProjectiveRotation(prefix)); + motherVolume->AddNode(volume, 1, transform); + }; + + // ------------------------------------------------------------ + // Optional trapezoidal frame + // ------------------------------------------------------------ + // This is exactly the old photosensor envelope. It is created for + // reference, but deliberately not added to the geometry. + photoFrames[photTileCount] = makeFrameFootprint(photThick); + auto* photoFrameVol = new TGeoVolume(Form("photoFrame_%d_%d", rPosId, photTileCount), photoFrames[photTileCount], medSi); + photoFrameVol->SetLineColor(kGray + 2); + photoFrameVol->SetLineWidth(1); + // Uncomment only when the mechanical frame material/solid geometry should be included. + // This would be a solid trapezoid and would overlap the sensitive silicon: need for opening + // motherVolume->AddNode(photoFrameVol, 1, new TGeoCombiTrans(photoCenterR * TMath::Cos(phiRad), photoCenterR * TMath::Sin(phiRad), photoCenterZ, makeProjectiveRotation("photoFrame"))); + + // ------------------------------------------------------------ + // True sensitive silicon: centered 17 x 18 cm2 rectangle + // ------------------------------------------------------------ + // Local X corresponds to the in-plane Z direction after placement. + // Local Y corresponds to the in-plane r-phi direction. + // Local Z is the 1 mm thickness direction. + /*photoTile = new TGeoArb8(photThick / 2.0); + photoTile->SetVertex(0, sipmActiveSizeZ / 2.0, -sipmActiveSizeRPhi / 2.0); + photoTile->SetVertex(1, -sipmActiveSizeZ / 2.0, -sipmActiveSizeRPhi / 2.0); + photoTile->SetVertex(2, -sipmActiveSizeZ / 2.0, sipmActiveSizeRPhi / 2.0); + photoTile->SetVertex(3, sipmActiveSizeZ / 2.0, sipmActiveSizeRPhi / 2.0); + photoTile->SetVertex(4, sipmActiveSizeZ / 2.0, -sipmActiveSizeRPhi / 2.0); + photoTile->SetVertex(5, -sipmActiveSizeZ / 2.0, -sipmActiveSizeRPhi / 2.0); + photoTile->SetVertex(6, -sipmActiveSizeZ / 2.0, sipmActiveSizeRPhi / 2.0); + photoTile->SetVertex(7, sipmActiveSizeZ / 2.0, sipmActiveSizeRPhi / 2.0); + auto* photoTileVol = new TGeoVolume(Form("%s_%d_%d", GeometryTGeo::getRICHSensorPattern(), rPosId, photTileCount), photoTile, medSi); + photoTileVol->SetLineColor(kRed); + photoTileVol->SetLineWidth(1); + auto* rotTransPhoto = new TGeoCombiTrans(photoCenterR * TMath::Cos(phiRad), photoCenterR * TMath::Sin(phiRad), photoCenterZ, makeProjectiveRotation("photoTile")); + motherVolume->AddNode(photoTileVol, 1, rotTransPhoto);*/ + + // Silicone resin in front of SiPMs + addReadoutLayer("siliconeLayer", siliconeLayerThickness, siliconeCenterOffset, medSilicone, kOrange + 2, useRectangularModules); + + // Active sensitive silicon. + photoTile = makeRectangularFootprint(activeSiliconThickness); + auto* photoTileVol = new TGeoVolume(Form("%s_%d_%d", GeometryTGeo::getRICHSensorPattern(), rPosId, photTileCount), photoTile, medSi); + const double activeSiliconCenterR = photoCenterR + activeSiliconCenterOffset * normalRadial; + const double activeSiliconCenterZ = photoCenterZ + activeSiliconCenterOffset * normalZ; + auto* rotTransPhoto = new TGeoCombiTrans(activeSiliconCenterR * TMath::Cos(phiRad), activeSiliconCenterR * TMath::Sin(phiRad), activeSiliconCenterZ, makeProjectiveRotation("photoTile")); + motherVolume->AddNode(photoTileVol, 1, rotTransPhoto); + + // Passive silicon absorber. + addReadoutLayer("siliconAbsorber", passiveSiliconThickness, passiveSiliconCenterOffset, medSiAbsorber, kBlue + 1, true); + + // ------------------------------------------------------------ + // Stack behind the SiPM + // ------------------------------------------------------------ + // Every gap is surface-to-surface. centerOffset is measured from + // the SiPM center along the outward local normal. + double outerSurfaceOffset = photThick / 2.0; + + outerSurfaceOffset += gapSiPMToPCB1; + const double pcb1CenterOffset = outerSurfaceOffset + pcb1Thickness / 2.0; + addReadoutLayer("pcb1", pcb1Thickness, pcb1CenterOffset, medFR4, kGreen + 1, useRectangularModules); + outerSurfaceOffset += pcb1Thickness; + + outerSurfaceOffset += gapPCB1ToCoolingPlate; + const double coolingPlateCenterOffset = outerSurfaceOffset + coolingPlateThickness / 2.0; + addReadoutLayer("coolingPlate", coolingPlateThickness, coolingPlateCenterOffset, medHTCC, kRed, useRectangularModules); + outerSurfaceOffset += coolingPlateThickness; + + outerSurfaceOffset += gapCoolingPlateToPCB2; + const double pcb2CenterOffset = outerSurfaceOffset + pcb2Thickness / 2.0; + addReadoutLayer("pcb2", pcb2Thickness, pcb2CenterOffset, medFR4, kGreen + 2, useRectangularModules); + outerSurfaceOffset += pcb2Thickness; + + outerSurfaceOffset += gapPCB2ToPCB3; + const double pcb3CenterOffset = outerSurfaceOffset + pcb3Thickness / 2.0; + addReadoutLayer("pcb3", pcb3Thickness, pcb3CenterOffset, medFR4, kGreen + 3, useRectangularModules); + + photTileCount++; + } + } + + // Gas sectors (argon) - legacy code, not used in the current geometry, but kept for reference + /* + for (auto& gasSector : gasSectors) { + double separation{(aerDetDistance - radThick - photThick)}; auto* radiator = radiatorTiles[argSectorsCount]; auto* photosensor = photoTiles[argSectorsCount]; - argonSector = new TGeoArb8(separation / 2); - - argonSector->SetVertex(0, -photZ / 2, -photYmin / 2); - argonSector->SetVertex(1, -photZ / 2, photYmin / 2); - argonSector->SetVertex(2, photZ / 2, photYmax / 2); - argonSector->SetVertex(3, photZ / 2, -photYmax / 2); - argonSector->SetVertex(4, -radZ / 2, -radYmin / 2); - argonSector->SetVertex(5, -radZ / 2, radYmin / 2); - argonSector->SetVertex(6, radZ / 2, radYmax / 2); - argonSector->SetVertex(7, radZ / 2, -radYmax / 2); - - TGeoVolume* argonSectorVol = new TGeoVolume(Form("argonSector_%d_%d", rPosId, argSectorsCount), argonSector, medAr); - argonSectorVol->SetVisibility(kTRUE); - argonSectorVol->SetLineColor(kOrange - 8); - argonSectorVol->SetLineWidth(1); - auto* rotArgon = new TGeoRotation(Form("argonSectorRotation_%d_%d", argSectorsCount, rPosId)); - rotArgon->RotateY(-90 - thetaBDeg); - rotArgon->RotateZ(argSectorsCount * deltaPhiDeg); - auto* rotTransArgon = new TGeoCombiTrans((radRad0 + TMath::Cos(thetaB) * (separation + radThick) / 2) * TMath::Cos(argSectorsCount * TMath::Pi() / (nTilesPhi / 2)), - (radRad0 + TMath::Cos(thetaB) * (separation + radThick) / 2) * TMath::Sin(argSectorsCount * TMath::Pi() / (nTilesPhi / 2)), - radRad0 * TMath::Tan(thetaB) + TMath::Sin(thetaB) * (separation + radThick) / 2, - rotArgon); - motherVolume->AddNode(argonSectorVol, 1, rotTransArgon); + gasSector = new TGeoArb8(separation / 2); + + gasSector->SetVertex(0, -photZ / 2, -photYmin / 2); + gasSector->SetVertex(1, -photZ / 2, photYmin / 2); + gasSector->SetVertex(2, photZ / 2, photYmax / 2); + gasSector->SetVertex(3, photZ / 2, -photYmax / 2); + gasSector->SetVertex(4, -radZ / 2, -radYmin / 2); + gasSector->SetVertex(5, -radZ / 2, radYmin / 2); + gasSector->SetVertex(6, radZ / 2, radYmax / 2); + gasSector->SetVertex(7, radZ / 2, -radYmax / 2); + + TGeoVolume* gasSectorVol = new TGeoVolume(Form("gasSector_%d_%d", rPosId, argSectorsCount), gasSector, medCO2); + gasSectorVol->SetVisibility(kTRUE); + gasSectorVol->SetLineColor(kOrange - 8); + gasSectorVol->SetLineWidth(1); + auto* rotGas = new TGeoRotation(Form("gasSectorRotation_%d_%d", argSectorsCount, rPosId)); + rotGas->RotateY(-90 - thetaBDeg); + //rotGas->RotateZ(argSectorsCount * deltaPhiDeg); + //auto* rotTransGas = new TGeoCombiTrans((radRad0 + TMath::Cos(thetaB) * (separation + radThick) / 2) * TMath::Cos(argSectorsCount * TMath::Pi() / (nTilesPhi / 2)), + // (radRad0 + TMath::Cos(thetaB) * (separation + radThick) / 2) * TMath::Sin(argSectorsCount * TMath::Pi() / (nTilesPhi / 2)), + // radRad0 * TMath::Tan(thetaB) + TMath::Sin(thetaB) * (separation + radThick) / 2, + // rotGas); + const double gasPhiRad = modulePhiRad(argSectorsCount); + rotGas->RotateZ(gasPhiRad * 180.0 / TMath::Pi()); + auto* rotTransGas = new TGeoCombiTrans((radRad0 + TMath::Cos(thetaB) * (separation + radThick) / 2.0) * TMath::Cos(gasPhiRad), + (radRad0 + TMath::Cos(thetaB) * (separation + radThick) / 2.0) * TMath::Sin(gasPhiRad), + radRad0 * TMath::Tan(thetaB) + TMath::Sin(thetaB) * (separation + radThick) / 2.0, + rotGas); + motherVolume->AddNode(gasSectorVol, 1, rotTransGas); argSectorsCount++; } + */ } FWDRich::FWDRich(std::string name, - float rMin, - float rMax, - float zAerogelMin, - float dZAerogel, - float zArgonMin, - float dZArgon, - float zSiliconMin, - float dZSilicon) : mName{name}, - mRmin{rMin}, - mRmax{rMax}, - mZAerogelMin{zAerogelMin}, - mDZAerogel{dZAerogel}, - mZArgonMin{zArgonMin}, - mDZArgon{dZArgon}, - mZSiliconMin{zSiliconMin}, - mDZSilicon{dZSilicon} + double rMin, + double rMax, + double zAerogelMin, + double dZAerogel, + double zArgonMin, + double dZArgon, + double zSiliconMin, + double dZSilicon) : mName{name}, + mRmin{rMin}, + mRmax{rMax}, + mZAerogelMin{zAerogelMin}, + mDZAerogel{dZAerogel}, + mZArgonMin{zArgonMin}, + mDZArgon{dZArgon}, + mZSiliconMin{zSiliconMin}, + mDZSilicon{dZSilicon} { } BWDRich::BWDRich(std::string name, - float rMin, - float rMax, - float zAerogelMin, - float dZAerogel, - float zArgonMin, - float dZArgon, - float zSiliconMin, - float dZSilicon) : mName{name}, - mRmin{rMin}, - mRmax{rMax}, - mZAerogelMin{zAerogelMin}, - mDZAerogel{dZAerogel}, - mZArgonMin{zArgonMin}, - mDZArgon{dZArgon}, - mZSiliconMin{zSiliconMin}, - mDZSilicon{dZSilicon} + double rMin, + double rMax, + double zAerogelMin, + double dZAerogel, + double zArgonMin, + double dZArgon, + double zSiliconMin, + double dZSilicon) : mName{name}, + mRmin{rMin}, + mRmax{rMax}, + mZAerogelMin{zAerogelMin}, + mDZAerogel{dZAerogel}, + mZArgonMin{zArgonMin}, + mDZArgon{dZArgon}, + mZSiliconMin{zSiliconMin}, + mDZSilicon{dZSilicon} { } From 5e3551a597129702e281ff59a8687eeb0b948bc2 Mon Sep 17 00:00:00 2001 From: Nicola Nicassio Date: Fri, 17 Jul 2026 20:08:15 +0200 Subject: [PATCH 2/4] Latest v3b.1 RICH geometry with quadrants, modules and shielding --- .../ALICE3/RICH/base/include/RICHBase/RICHBaseParam.h | 10 +++++----- .../Upgrades/ALICE3/RICH/simulation/src/Detector.cxx | 2 +- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/Detectors/Upgrades/ALICE3/RICH/base/include/RICHBase/RICHBaseParam.h b/Detectors/Upgrades/ALICE3/RICH/base/include/RICHBase/RICHBaseParam.h index 52f5725ed032a..d53dab7756c8e 100644 --- a/Detectors/Upgrades/ALICE3/RICH/base/include/RICHBase/RICHBaseParam.h +++ b/Detectors/Upgrades/ALICE3/RICH/base/include/RICHBase/RICHBaseParam.h @@ -28,20 +28,20 @@ struct RICHBaseParam : public o2::conf::ConfigurableParamHelper { int nRings = 11; // (25 in v3) int nTiles = 44; // (36 in v3) bool oddGeom = true; // (false in v3) - + // The active and passive silicon thicknesses must sum to detectorThickness. double siliconeLayerThickness = 0.010; // cm: 0.1 mm resin layer in front double detectorThickness = 0.1; // cm double activeSiliconThickness = 0.01; // cm: 0.1 mm sensitive silicon //double passiveSiliconThickness = 0.09f; // cm: (detectorThickness - activeSiliconThickness) - + //cylindrical aerogel layout bool useCylindricalAerogel = true; double cylindricalAerogelEtaRef = 0.85; - + // Enable geometry with rectangular modules bool useRectangularModules = true; - + // Barrel photosensor active area. double sipmActiveSizeZ = 18.0; // cm double sipmActiveSizeRPhi = 17.0; // cm @@ -71,7 +71,7 @@ struct RICHBaseParam : public o2::conf::ConfigurableParamHelper { double gapPCB1ToCoolingPlate = 0.10; // cm double gapCoolingPlateToPCB2 = 0.10; // cm double gapPCB2ToPCB3 = 0.10; // cm - + // Minimum edge-to-edge clearances used to avoid exact contacts between adjacent modules. double moduleClearanceZ = 0.02; // cm double moduleClearanceRPhi = 0.02; // cm diff --git a/Detectors/Upgrades/ALICE3/RICH/simulation/src/Detector.cxx b/Detectors/Upgrades/ALICE3/RICH/simulation/src/Detector.cxx index fd415b76564de..0dd938b931665 100644 --- a/Detectors/Upgrades/ALICE3/RICH/simulation/src/Detector.cxx +++ b/Detectors/Upgrades/ALICE3/RICH/simulation/src/Detector.cxx @@ -619,7 +619,7 @@ void Detector::createGeometry() outerWallQuadrantVolume->SetLineColor(kOrange - 8); // kGray endCapPlusQuadrantVolume->SetLineColor(kOrange - 8); // kGray endCapMinusQuadrantVolume->SetLineColor(kOrange - 8); // kGray - lateralStartWallVolume->SetLineColor(kOrange - 8); // kGray + lateralStartWallVolume->SetLineColor(kOrange - 8); // kGray lateralEndWallVolume->SetLineColor(kOrange - 8); // kGray innerWallQuadrantVolume->SetTransparency(0); // 80 From 1339c9b5bad59c818ce6fe2d4a7744d9d4777b16 Mon Sep 17 00:00:00 2001 From: Nicola Nicassio Date: Fri, 17 Jul 2026 20:18:32 +0200 Subject: [PATCH 3/4] Latest v3b.1 RICH geometry with quadrants, modules and shielding (clang) --- .../RICH/base/include/RICHBase/RICHBaseParam.h | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/Detectors/Upgrades/ALICE3/RICH/base/include/RICHBase/RICHBaseParam.h b/Detectors/Upgrades/ALICE3/RICH/base/include/RICHBase/RICHBaseParam.h index d53dab7756c8e..2d1d83d376ea3 100644 --- a/Detectors/Upgrades/ALICE3/RICH/base/include/RICHBase/RICHBaseParam.h +++ b/Detectors/Upgrades/ALICE3/RICH/base/include/RICHBase/RICHBaseParam.h @@ -25,17 +25,17 @@ struct RICHBaseParam : public o2::conf::ConfigurableParamHelper { double rMin = 104.0; // cm (90.0 in v3) double radiatorThickness = 2.0; // cm double zRichLength = 700.0; // cm - int nRings = 11; // (25 in v3) - int nTiles = 44; // (36 in v3) - bool oddGeom = true; // (false in v3) + int nRings = 11; // (25 in v3) + int nTiles = 44; // (36 in v3) + bool oddGeom = true; // (false in v3) // The active and passive silicon thicknesses must sum to detectorThickness. - double siliconeLayerThickness = 0.010; // cm: 0.1 mm resin layer in front - double detectorThickness = 0.1; // cm + double siliconeLayerThickness = 0.010; // cm: 0.1 mm resin layer in front + double detectorThickness = 0.1; // cm double activeSiliconThickness = 0.01; // cm: 0.1 mm sensitive silicon - //double passiveSiliconThickness = 0.09f; // cm: (detectorThickness - activeSiliconThickness) + // double passiveSiliconThickness = 0.09f; // cm: (detectorThickness - activeSiliconThickness) - //cylindrical aerogel layout + // cylindrical aerogel layout bool useCylindricalAerogel = true; double cylindricalAerogelEtaRef = 0.85; @@ -110,4 +110,4 @@ struct RICHBaseParam : public o2::conf::ConfigurableParamHelper { } // namespace rich } // end namespace o2 -#endif \ No newline at end of file +#endif From 612b9158edbbb3fce1f1bdacbec52b431642b4c6 Mon Sep 17 00:00:00 2001 From: Nicola Nicassio Date: Sat, 18 Jul 2026 01:23:49 +0200 Subject: [PATCH 4/4] Latest v3b.1 RICH geometry with quadrants, modules and shielding (clang fixed) --- .../RICH/simulation/include/RICHSimulation/RICHRing.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Detectors/Upgrades/ALICE3/RICH/simulation/include/RICHSimulation/RICHRing.h b/Detectors/Upgrades/ALICE3/RICH/simulation/include/RICHSimulation/RICHRing.h index 18c4c4e2f8d71..a1b8974c0b748 100644 --- a/Detectors/Upgrades/ALICE3/RICH/simulation/include/RICHSimulation/RICHRing.h +++ b/Detectors/Upgrades/ALICE3/RICH/simulation/include/RICHSimulation/RICHRing.h @@ -60,9 +60,9 @@ class Ring private: int mPosId; // id of the ring int mNTiles; // number of modules - double mRRad; // max distance for radiators - double mRPhot; // max distance for photosensitive surfaces - double mRadThickness; // thickness of the radiator + double mRRad; // max distance for radiators + double mRPhot; // max distance for photosensitive surfaces + double mRadThickness; // thickness of the radiator double mPhotThickness; // thickness of the photosensitive surface ClassDef(Ring, 0);