Back to home page

EIC code displayed by LXR
 
 

    

Warning, file /epic/src/DRICH_geo.cpp was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 - 2026 Christopher Dilks, Junhuai Xu, Luisa Occhiuto
0003 
0004 // - Design Adapted from Standalone Fun4all and GEMC implementations
0005 //   [ Evaristo Cisbani, Cristiano Fanelli, Alessio Del Dotto, et al. ]
0006 
0007 #include "DD4hep/DetFactoryHelper.h"
0008 #include "DD4hep/OpticalSurfaces.h"
0009 #include "DD4hep/Printout.h"
0010 #include "DDRec/DetectorData.h"
0011 #include "DDRec/Surface.h"
0012 #include <numeric>
0013 #include <vector>
0014 
0015 #include <XML/Helper.h>
0016 
0017 using namespace dd4hep;
0018 using namespace dd4hep::rec;
0019 
0020 #ifdef WITH_IRT2_SUPPORT
0021 #include <TFile.h>
0022 
0023 #include "IRT2/CherenkovDetectorCollection.h"
0024 #include "IRT2/SphericalSurface.h"
0025 
0026 using namespace IRT2;
0027 #endif
0028 
0029 // create the detector
0030 static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetector sens) {
0031 
0032   xml::DetElement detElem       = handle;
0033   std::string detName           = detElem.nameStr();
0034   int detID                     = detElem.id();
0035   xml::Component dims           = detElem.dimensions();
0036   OpticalSurfaceManager surfMgr = desc.surfaceManager();
0037   DetElement det(detName, detID);
0038   sens.setType("tracker");
0039 
0040   // Start optical configuration if needed;
0041 #ifdef WITH_IRT2_SUPPORT
0042   auto geometry = CherenkovDetectorCollection::Instance();
0043   auto cdet     = geometry->AddNewDetector(detName.c_str());
0044 #endif
0045 
0046   // attributes, from compact file =============================================
0047   // - vessel
0048   auto vesselZmin   = dims.attr<double>(_Unicode(zmin));
0049   auto vesselLength = dims.attr<double>(_Unicode(length));
0050   auto vesselRmin0  = dims.attr<double>(_Unicode(rmin0));
0051   auto vesselRmin1  = dims.attr<double>(_Unicode(rmin1));
0052   auto vesselRmax0  = dims.attr<double>(_Unicode(rmax0));
0053   auto vesselRmax1  = dims.attr<double>(_Unicode(rmax1));
0054   auto vesselRmax2  = dims.attr<double>(_Unicode(rmax2));
0055   auto snoutLength  = dims.attr<double>(_Unicode(snout_length));
0056   auto nSectors     = dims.attr<int>(_Unicode(nsectors));
0057 #ifdef WITH_IRT2_SUPPORT
0058   cdet->SetSectorCount(nSectors);
0059   // The way Chris defined it in the geometry;
0060   cdet->SetSectorPhase(0.0);
0061 #endif
0062   auto wallThickness   = dims.attr<double>(_Unicode(wall_thickness));
0063   auto windowThickness = dims.attr<double>(_Unicode(window_thickness));
0064   auto vesselMat       = desc.material(detElem.attr<std::string>(_Unicode(material)));
0065 
0066   auto gasvolMatName = detElem.attr<std::string>(_Unicode(gas));
0067   auto gasvolMat     = desc.material(gasvolMatName);
0068 
0069   auto vesselVis = desc.visAttributes(detElem.attr<std::string>(_Unicode(vis_vessel)));
0070   auto gasvolVis = desc.visAttributes(detElem.attr<std::string>(_Unicode(vis_gas)));
0071   // - radiator (applies to aerogel and filter)
0072   auto radiatorElem       = detElem.child(_Unicode(radiator));
0073   auto radiatorRmin       = radiatorElem.attr<double>(_Unicode(rmin));
0074   auto radiatorRmax       = radiatorElem.attr<double>(_Unicode(rmax));
0075   auto radiatorPitch      = radiatorElem.attr<double>(_Unicode(pitch));
0076   auto radiatorFrontplane = radiatorElem.attr<double>(_Unicode(frontplane));
0077   // - aerogel
0078   auto aerogelElem      = radiatorElem.child(_Unicode(aerogel));
0079   auto aerogelMatName   = aerogelElem.attr<std::string>(_Unicode(material));
0080   auto aerogelMat       = desc.material(aerogelMatName);
0081   auto aerogelVis       = desc.visAttributes(aerogelElem.attr<std::string>(_Unicode(vis)));
0082   auto aerogelThickness = aerogelElem.attr<double>(_Unicode(thickness));
0083   // - aerogel structure
0084   auto coronasElem      = radiatorElem.child(_Unicode(coronas));
0085   auto coronasMat       = desc.material(coronasElem.attr<std::string>(_Unicode(material)));
0086   auto coronasVis       = desc.visAttributes(coronasElem.attr<std::string>(_Unicode(vis)));
0087   auto coronasThickness = coronasElem.attr<double>(_Unicode(thickness));
0088   auto segmentationType = coronasElem.attr<std::string>("segmentation");
0089   // read crown parameters from child <crown> elements
0090   std::vector<double> radii;
0091   std::vector<int> numSegments;
0092   for (xml::Collection_t crownIt(coronasElem, _Unicode(crown)); crownIt; ++crownIt) {
0093     xml::Component crownElem = crownIt;
0094     radii.push_back(crownElem.attr<double>(_Unicode(radius)));
0095     if (crownElem.hasAttr(_Unicode(num_segments)))
0096       numSegments.push_back(crownElem.attr<int>(_Unicode(num_segments)));
0097   }
0098   int numCrowns = radii.size();
0099   // - filter
0100   auto filterElem      = radiatorElem.child(_Unicode(filter));
0101   auto filterMatName   = filterElem.attr<std::string>(_Unicode(material));
0102   auto filterMat       = desc.material(filterMatName);
0103   auto filterVis       = desc.visAttributes(filterElem.attr<std::string>(_Unicode(vis)));
0104   auto filterThickness = filterElem.attr<double>(_Unicode(thickness));
0105   // - airgap between filter and aerogel
0106   auto airgapElem      = radiatorElem.child(_Unicode(airgap));
0107   auto airgapMat       = desc.material(airgapElem.attr<std::string>(_Unicode(material)));
0108   auto airgapVis       = desc.visAttributes(airgapElem.attr<std::string>(_Unicode(vis)));
0109   auto airgapThickness = airgapElem.attr<double>(_Unicode(thickness));
0110   // - mirror
0111   auto mirrorElem      = detElem.child(_Unicode(mirror));
0112   auto mirrorMat       = desc.material(mirrorElem.attr<std::string>(_Unicode(material)));
0113   auto mirrorVis       = desc.visAttributes(mirrorElem.attr<std::string>(_Unicode(vis)));
0114   auto mirrorSurf      = surfMgr.opticalSurface(mirrorElem.attr<std::string>(_Unicode(surface)));
0115   auto mirrorBackplane = mirrorElem.attr<double>(_Unicode(backplane));
0116   auto mirrorThickness = mirrorElem.attr<double>(_Unicode(thickness));
0117   auto mirrorRmin      = mirrorElem.attr<double>(_Unicode(rmin));
0118   auto mirrorRmax      = mirrorElem.attr<double>(_Unicode(rmax));
0119   auto mirrorPhiw      = mirrorElem.attr<double>(_Unicode(phiw));
0120   auto focusTuneZ      = mirrorElem.attr<double>(_Unicode(focus_tune_z));
0121   auto focusTuneX      = mirrorElem.attr<double>(_Unicode(focus_tune_x));
0122   // - sensorboxes
0123   auto sensorboxLength = desc.constant<double>("DRICH_sensorbox_length");
0124   auto sensorboxRmin   = desc.constant<double>("DRICH_sensorbox_rmin");
0125   auto sensorboxRmax   = desc.constant<double>("DRICH_sensorbox_rmax");
0126   auto sensorboxDphi   = desc.constant<double>("DRICH_sensorbox_dphi");
0127   // - sensor photosensitive surface (pss)
0128   auto pssElem      = detElem.child(_Unicode(sensors)).child(_Unicode(pss));
0129   auto pssMat       = desc.material(pssElem.attr<std::string>(_Unicode(material)));
0130   auto pssVis       = desc.visAttributes(pssElem.attr<std::string>(_Unicode(vis)));
0131   auto pssSurf      = surfMgr.opticalSurface(pssElem.attr<std::string>(_Unicode(surface)));
0132   auto pssSide      = pssElem.attr<double>(_Unicode(side));
0133   auto pssThickness = pssElem.attr<double>(_Unicode(thickness));
0134   // - sensor resin
0135   auto resinElem      = detElem.child(_Unicode(sensors)).child(_Unicode(resin));
0136   auto resinMat       = desc.material(resinElem.attr<std::string>(_Unicode(material)));
0137   auto resinVis       = desc.visAttributes(resinElem.attr<std::string>(_Unicode(vis)));
0138   auto resinSide      = resinElem.attr<double>(_Unicode(side));
0139   auto resinThickness = resinElem.attr<double>(_Unicode(thickness));
0140   // - photodetector unit (PDU)
0141   auto pduElem       = detElem.child(_Unicode(sensors)).child(_Unicode(pdu));
0142   auto pduNumSensors = desc.constant<int>("DRICH_pdu_num_sensors");
0143   auto pduSensorGap  = desc.constant<double>("DRICH_pdu_sensor_gap");
0144   auto pduGap        = desc.constant<double>("DRICH_pdu_gap");
0145   // - sensor sphere
0146   auto sensorSphElem    = detElem.child(_Unicode(sensors)).child(_Unicode(sphere));
0147   auto sensorSphRadius  = sensorSphElem.attr<double>(_Unicode(radius));
0148   auto sensorSphCenterX = sensorSphElem.attr<double>(_Unicode(centerx));
0149   auto sensorSphCenterZ = sensorSphElem.attr<double>(_Unicode(centerz));
0150   // - sensor sphere patch cuts
0151   auto sensorSphPatchElem = detElem.child(_Unicode(sensors)).child(_Unicode(sphericalpatch));
0152   auto sensorSphPatchPhiw = sensorSphPatchElem.attr<double>(_Unicode(phiw));
0153   auto sensorSphPatchRmin = sensorSphPatchElem.attr<double>(_Unicode(rmin));
0154   auto sensorSphPatchRmax = sensorSphPatchElem.attr<double>(_Unicode(rmax));
0155   auto sensorSphPatchZmin = sensorSphPatchElem.attr<double>(_Unicode(zmin));
0156   // - sensor readout
0157   auto readoutName = detElem.attr<std::string>(_Unicode(readout));
0158   // - settings and switches
0159   auto debugOpticsMode = desc.constant<int>("DRICH_debug_optics");
0160   bool debugSector     = desc.constant<int>("DRICH_debug_sector") == 1;
0161   bool debugMirror     = desc.constant<int>("DRICH_debug_mirror") == 1;
0162   bool debugSensors    = desc.constant<int>("DRICH_debug_sensors") == 1;
0163 
0164   // if debugging optics, override some settings
0165   bool debugOptics = debugOpticsMode > 0;
0166   if (debugOptics) {
0167     printout(WARNING, "DRICH_geo", "DEBUGGING DRICH OPTICS");
0168     switch (debugOpticsMode) {
0169     case 1:
0170       vesselMat = aerogelMat = filterMat = pssMat = gasvolMat = desc.material("VacuumOptical");
0171       break;
0172     case 2:
0173       vesselMat = aerogelMat = filterMat = pssMat = desc.material("VacuumOptical");
0174       break;
0175     case 3:
0176       vesselMat = aerogelMat = filterMat = gasvolMat = desc.material("VacuumOptical");
0177       break;
0178     default:
0179       printout(FATAL, "DRICH_geo", "UNKNOWN debugOpticsMode");
0180       return det;
0181     }
0182   }
0183 
0184   // if debugging anything, draw only one sector and adjust visibility
0185   if (debugOptics || debugMirror || debugSensors)
0186     debugSector = true;
0187   if (debugSector)
0188     gasvolVis = vesselVis = desc.invisible();
0189 
0190 #ifdef WITH_IRT2_SUPPORT
0191   SphericalSurface* msurface       = 0;
0192   IRT2::OpticalBoundary* mboundary = 0;
0193 #endif
0194 
0195   // readout coder <-> unique sensor ID
0196   /* - `sensorIDfields` is a list of readout fields used to specify a unique sensor ID
0197    * - `cellMask` is defined such that a hit's `cellID & cellMask` is the corresponding sensor's unique ID
0198    */
0199   std::vector<std::string> sensorIDfields = {"pdu", "sipm", "sector"};
0200   const auto& readoutCoder                = *desc.readout(readoutName).idSpec().decoder();
0201   // determine `cellMask` based on `sensorIDfields`
0202   uint64_t cellMask = 0;
0203   for (const auto& idField : sensorIDfields)
0204     cellMask |= readoutCoder[idField].mask();
0205   desc.add(Constant("DRICH_cell_mask", std::to_string(cellMask)));
0206   // create a unique sensor ID from a sensor's PlacedVolume::volIDs
0207   auto encodeSensorID = [&readoutCoder](auto ids) {
0208     uint64_t enc = 0;
0209     for (const auto& [idField, idValue] : ids)
0210       enc |= uint64_t(idValue) << readoutCoder[idField].offset();
0211     return enc;
0212   };
0213 #ifdef WITH_IRT2_SUPPORT
0214   uint64_t sector_mask = ~(0x0ul) ^ readoutCoder["sector"].mask();
0215   // Want to mask away dRICH sector bits in this mask;
0216   cdet->SetReadoutCellMask(cellMask & sector_mask);
0217 #endif
0218 
0219   // BUILD VESSEL ====================================================================
0220   /* - `vessel`: aluminum enclosure, the mother volume of the dRICH
0221    * - `gasvol`: gas volume, which fills `vessel`; all other volumes defined below
0222    *   are children of `gasvol`
0223    * - the dRICH vessel geometry has two regions: the snout refers to the conic region
0224    *   in the front, housing the aerogel, while the tank refers to the cylindrical
0225    *   region, housing the rest of the detector components
0226    */
0227 
0228   // derived attributes
0229   double tankLength = vesselLength - snoutLength;
0230   double vesselZmax = vesselZmin + vesselLength;
0231 
0232   // snout solids
0233   double boreDelta  = vesselRmin1 - vesselRmin0;
0234   double snoutDelta = vesselRmax1 - vesselRmax0;
0235   Cone vesselSnout(snoutLength / 2.0, vesselRmin0, vesselRmax0,
0236                    vesselRmin0 + boreDelta * snoutLength / vesselLength, vesselRmax1);
0237   Cone gasvolSnout(
0238       /* note: `gasvolSnout` extends a bit into the tank, so it touches `gasvolTank`
0239        * - the extension distance is equal to the tank `windowThickness`, so the
0240        *   length of `gasvolSnout` == length of `vesselSnout`
0241        * - the extension backplane radius is calculated using similar triangles
0242        */
0243       snoutLength / 2.0, vesselRmin0 + wallThickness, vesselRmax0 - wallThickness,
0244       vesselRmin0 + boreDelta * (snoutLength - windowThickness) / vesselLength + wallThickness,
0245       vesselRmax1 - wallThickness + windowThickness * (vesselRmax1 - vesselRmax0) / snoutLength);
0246 
0247   // tank solids:
0248   // - inner: cone along beamline
0249   // - outer: cone to back of sensor box, then fixed radius cylinder
0250   Polycone vesselTank(
0251       0, 2 * M_PI,
0252       /* rmin */
0253       {vesselSnout.rMin2(),
0254        std::lerp(vesselSnout.rMin2(), vesselRmin1, (sensorboxLength - snoutLength) / tankLength),
0255        vesselRmin1},
0256       /* rmax */ {vesselSnout.rMax2(), vesselRmax2, vesselRmax2},
0257       /* z    */
0258       {-tankLength / 2.0, -tankLength / 2.0 + sensorboxLength - snoutLength, tankLength / 2.0});
0259   Polycone gasvolTank(
0260       0, 2 * M_PI,
0261       /* rmin */
0262       {gasvolSnout.rMin2(),
0263        std::lerp(gasvolSnout.rMin2(), vesselRmin1 + wallThickness,
0264                  (sensorboxLength - snoutLength) / tankLength),
0265        vesselRmin1 + wallThickness},
0266       /* rmax */ {gasvolSnout.rMax2(), vesselRmax2 - wallThickness, vesselRmax2 - wallThickness},
0267       /* z    */
0268       {-tankLength / 2.0 + windowThickness,
0269        -tankLength / 2.0 + windowThickness + sensorboxLength - snoutLength,
0270        tankLength / 2.0 - windowThickness});
0271 
0272   // sensorbox solids
0273   double dphi = atan2(wallThickness, sensorboxRmax); // thickness only correct at Rmax
0274   Tube vesselSensorboxTube(sensorboxRmin, sensorboxRmax, sensorboxLength / 2., -sensorboxDphi / 2.,
0275                            sensorboxDphi / 2.);
0276   Tube gasvolSensorboxTube(sensorboxRmin + wallThickness, sensorboxRmax - wallThickness,
0277                            sensorboxLength / 2., -sensorboxDphi / 2. + dphi,
0278                            sensorboxDphi / 2. - dphi);
0279 
0280   // union: snout + tank
0281   UnionSolid vesselUnion(vesselTank, vesselSnout, Position(0., 0., -vesselLength / 2.));
0282   UnionSolid gasvolUnion(gasvolTank, gasvolSnout,
0283                          Position(0., 0., -vesselLength / 2. + windowThickness));
0284 
0285   // union: add sensorboxes for all sectors
0286   for (int isec = 0; isec < nSectors; isec++) {
0287     RotationZ sectorRotation((isec + 0.5) * 2 * M_PI / nSectors);
0288     vesselUnion = UnionSolid(
0289         vesselUnion, vesselSensorboxTube,
0290         Transform3D(sectorRotation, Position(0., 0., -(snoutLength + sensorboxLength - 0.6) / 2.)));
0291     gasvolUnion = UnionSolid(
0292         gasvolUnion, gasvolSensorboxTube,
0293         Transform3D(sectorRotation,
0294                     Position(0., 0., -(snoutLength + sensorboxLength) / 2. + windowThickness)));
0295   }
0296 
0297   //  extra solids for `debugOptics` only
0298   Box vesselBox(1001, 1001, 1001);
0299   Box gasvolBox(1000, 1000, 1000);
0300 
0301   // choose vessel and gasvol solids (depending on `debugOpticsMode` (0=disabled))
0302   Solid vesselSolid, gasvolSolid;
0303   switch (debugOpticsMode) {
0304   case 0:
0305     vesselSolid = vesselUnion;
0306     gasvolSolid = gasvolUnion;
0307     break; // `!debugOptics`
0308   case 1:
0309   case 3:
0310     vesselSolid = vesselBox;
0311     gasvolSolid = gasvolBox;
0312     break;
0313   case 2:
0314     vesselSolid = vesselBox;
0315     gasvolSolid = gasvolUnion;
0316     break;
0317   }
0318 
0319   // volumes
0320   Volume vesselVol(detName, vesselSolid, vesselMat);
0321   Volume gasvolVol(detName + "_gas", gasvolSolid, gasvolMat);
0322   vesselVol.setVisAttributes(vesselVis);
0323   vesselVol.setRegion(desc, detElem.regionStr());
0324   gasvolVol.setVisAttributes(gasvolVis);
0325 
0326   // reference positions
0327   // - the vessel is created such that the center of the cylindrical tank volume
0328   //   coincides with the origin; this is called the "origin position" of the vessel
0329   // - when the vessel (and its children volumes) is placed, it is translated in
0330   //   the z-direction to be in the proper EPIC-integration location
0331   // - these reference positions are for the frontplane and backplane of the vessel,
0332   //   with respect to the vessel origin position
0333   auto originFront = Position(0., 0., -tankLength / 2.0 - snoutLength);
0334   // auto originBack  = Position(0., 0., tankLength / 2.0);
0335   auto vesselPos = Position(0, 0, vesselZmin) - originFront;
0336 
0337   // place gas volume
0338   PlacedVolume gasvolPV = vesselVol.placeVolume(gasvolVol, Position(0, 0, 0));
0339   DetElement gasvolDE(det, "gasvol_de", 0);
0340   gasvolDE.setPlacement(gasvolPV);
0341 
0342   // place mother volume (vessel)
0343   Volume motherVol      = desc.pickMotherVolume(det);
0344   PlacedVolume vesselPV = motherVol.placeVolume(vesselVol, vesselPos);
0345   vesselPV.addPhysVolID("system", detID);
0346   det.setPlacement(vesselPV);
0347 
0348   // BUILD RADIATOR ====================================================================
0349 
0350   // solid and volume: create aerogel and filter
0351 
0352   // ------------------------------------------------------------------------
0353   Volume aerogelVol;
0354   PlacedVolume aerogelPV;
0355   DetElement aerogelDE;
0356   double structureThickness = aerogelThickness;
0357   auto radiatorPos = Position(0., 0., radiatorFrontplane + 0.5 * structureThickness) + originFront;
0358 
0359   Cone aerogelSolid(aerogelThickness / 2.0, radiatorRmin, radiatorRmax,
0360                     radiatorRmin + boreDelta * aerogelThickness / vesselLength, radiatorRmax);
0361   auto aerogelPlacement = Translation3D(radiatorPos) * // re-center to originFront
0362                           RotationY(radiatorPitch);    // change polar angle to specified pitch
0363   aerogelVol            = Volume(detName + "_aerogel", aerogelSolid, aerogelMat);
0364   aerogelVol.setVisAttributes(aerogelVis);
0365   aerogelPV = gasvolVol.placeVolume(aerogelVol, aerogelPlacement);
0366   aerogelDE = DetElement(det, "aerogel_de", 0);
0367   aerogelDE.setPlacement(aerogelPV);
0368 
0369   if (segmentationType == "trapezoidal") {
0370     double crownHeight = structureThickness;
0371     std::vector<double> innerRadiusBottoms_half;
0372     std::vector<double> innerRadiusTops_half;
0373     std::vector<double> outerRadiusBottoms_half;
0374     std::vector<double> outerRadiusTops_half;
0375 
0376     // Create and place individual crown volumes
0377     for (int i = 0; i < numCrowns; i++) {
0378       double centralRadius = radii[i];
0379       double rMinBottom, rMinTop, rMaxBottom, rMaxTop;
0380 
0381       if (i == 0) {
0382 
0383         rMinBottom = centralRadius - boreDelta * (coronasThickness / 2.0) / vesselLength;
0384         rMinTop    = rMinBottom + snoutDelta * aerogelThickness / snoutLength;
0385         rMaxBottom = rMinBottom + boreDelta * coronasThickness / vesselLength;
0386         rMaxTop    = rMinTop + boreDelta * coronasThickness / vesselLength;
0387 
0388       } else {
0389         double innerRadius = centralRadius - coronasThickness / 2.0;
0390         double outerRadius = centralRadius + coronasThickness / 2.0;
0391 
0392         // FIX PROTRUSION
0393 
0394         if (i == numCrowns - 1) {
0395           double safetyMargin = 0.05 * dd4hep::cm;
0396           outerRadius -= safetyMargin;
0397         }
0398 
0399         rMinBottom = innerRadius;
0400         rMinTop    = innerRadius;
0401         rMaxBottom = outerRadius;
0402         rMaxTop    = outerRadius;
0403       }
0404       innerRadiusBottoms_half.push_back(rMinBottom);
0405       innerRadiusTops_half.push_back(rMinTop);
0406       outerRadiusBottoms_half.push_back(rMaxBottom);
0407       outerRadiusTops_half.push_back(rMaxTop);
0408 
0409       Cone crownSolid(crownHeight / 2.0, rMinBottom, rMaxBottom, rMinTop, rMaxTop);
0410       std::string crownName = "CarbonCrown_" + std::to_string(i);
0411       Volume crownVol(crownName, crownSolid, coronasMat);
0412       crownVol.setVisAttributes(coronasVis);
0413 
0414       // Place crown volume directly in aerogel
0415       aerogelVol.placeVolume(crownVol, Position(0., 0., 0.));
0416     }
0417 
0418     // Create and place individual segment volumes
0419     for (int i = 0; i < numCrowns - 1; i++) {
0420       int N = numSegments[i];
0421 
0422       double rMin_Zminus = outerRadiusBottoms_half[i];
0423       double rMax_Zminus = innerRadiusBottoms_half[i + 1];
0424       double rMin_Zplus  = outerRadiusTops_half[i];
0425       double rMax_Zplus  = innerRadiusTops_half[i + 1];
0426 
0427       double segmentSpacing      = 2 * M_PI / N;
0428       double segmentAngularWidth = coronasThickness / rMin_Zminus;
0429 
0430       for (int p = 0; p < N; p++) {
0431         double phiStart = p * segmentSpacing;
0432         double phiEnd   = phiStart + segmentAngularWidth;
0433 
0434         ConeSegment segmentSolid(crownHeight / 2.0, rMin_Zminus, rMax_Zminus, rMin_Zplus,
0435                                  rMax_Zplus, phiStart, phiEnd);
0436         std::string segName = "CarbonSegment_" + std::to_string(i) + "_" + std::to_string(p);
0437         Volume segVol(segName, segmentSolid, coronasMat);
0438         segVol.setVisAttributes(coronasVis);
0439 
0440         // Place segment volume directly in aerogel
0441         aerogelVol.placeVolume(segVol, Position(0., 0., 0.));
0442       }
0443     } //crown
0444   } //trapezoidal
0445 
0446   else if (segmentationType == "square") {
0447     printout(WARNING, "DRICH_geo", "Square segmentation requested but not implemented yet.");
0448   } //Square
0449 
0450   Cone airgapSolid(airgapThickness / 2.0,
0451                    radiatorRmin + boreDelta * structureThickness / vesselLength, radiatorRmax,
0452                    radiatorRmin + boreDelta * (structureThickness + airgapThickness) / vesselLength,
0453                    radiatorRmax);
0454   Cone filterSolid(filterThickness / 2.0,
0455                    radiatorRmin + boreDelta * (structureThickness + airgapThickness) / vesselLength,
0456                    radiatorRmax,
0457                    radiatorRmin + boreDelta *
0458                                       (structureThickness + airgapThickness + filterThickness) /
0459                                       vesselLength,
0460                    radiatorRmax);
0461 
0462   Volume airgapVol(detName + "_airgap", airgapSolid, airgapMat);
0463   Volume filterVol(detName + "_filter", filterSolid, filterMat);
0464   airgapVol.setVisAttributes(airgapVis);
0465   filterVol.setVisAttributes(filterVis);
0466 
0467   // airgap and filter placement and surface properties
0468   if (!debugOptics) {
0469 
0470     auto airgapPlacement =
0471         Translation3D(radiatorPos) * // re-center to originFront
0472         RotationY(radiatorPitch) *   // change polar angle
0473         Translation3D(0., 0.,
0474                       (structureThickness + airgapThickness) / 2.); // move to aerogel backplane
0475     auto airgapPV = gasvolVol.placeVolume(airgapVol, airgapPlacement);
0476     DetElement airgapDE(det, "airgap_de", 0);
0477     airgapDE.setPlacement(airgapPV);
0478 
0479     auto filterPlacement =
0480         Translation3D(0., 0., airgapThickness) * // add an air gap
0481         Translation3D(radiatorPos) *             // re-center to originFront
0482         RotationY(radiatorPitch) *               // change polar angle
0483         Translation3D(0., 0.,
0484                       (structureThickness + filterThickness) / 2.); // move to aerogel backplane
0485     auto filterPV = gasvolVol.placeVolume(filterVol, filterPlacement);
0486     DetElement filterDE(det, "filter_de", 0);
0487     filterDE.setPlacement(filterPV);
0488 
0489 #if defined(WITH_IRT2_SUPPORT) || defined(WITH_IRT1_SUPPORT)
0490     // radiator z-positions (w.r.t. IP); only needed downstream if !debugOptics
0491     double aerogelZpos = vesselPos.z() + aerogelPV.position().z();
0492     double filterZpos  = vesselPos.z() + filterPV.position().z();
0493 #endif
0494 
0495 #ifdef WITH_IRT1_SUPPORT
0496     double airgapZpos = vesselPos.z() + airgapPV.position().z();
0497     desc.add(Constant("DRICH_aerogel_zpos", std::to_string(aerogelZpos)));
0498     desc.add(Constant("DRICH_airgap_zpos", std::to_string(airgapZpos)));
0499     desc.add(Constant("DRICH_filter_zpos", std::to_string(filterZpos)));
0500 #endif
0501 
0502 #ifdef WITH_IRT2_SUPPORT
0503     {
0504       TVector3 nx(1, 0, 0), ny(0, -1, 0);
0505 
0506       for (int isec = 0; isec < nSectors; isec++) {
0507         {
0508           // FIXME: Z-location does not really matter here, right?;
0509           auto boundary = new FlatSurface(TVector3(0, 0, 0), nx, ny);
0510 
0511           auto radiator = geometry->SetContainerVolume(cdet, "GasVolume", isec,
0512                                                        (G4LogicalVolume*)(0x0), 0, boundary);
0513           radiator->SetAlternativeMaterialName(gasvolMatName.c_str());
0514         }
0515 
0516         {
0517           auto surface = new FlatSurface((1 / mm) * TVector3(0, 0, aerogelZpos), nx, ny);
0518 
0519           auto radiator =
0520               geometry->AddFlatRadiator(cdet, "Aerogel", CherenkovDetector::Upstream, isec,
0521                                         (G4LogicalVolume*)(0x1), 0, surface, aerogelThickness / mm);
0522           radiator->SetAlternativeMaterialName(aerogelMatName.c_str());
0523         }
0524 
0525         {
0526           auto surface = new FlatSurface((1 / mm) * TVector3(0, 0, filterZpos), nx, ny);
0527 
0528           auto radiator =
0529               geometry->AddFlatRadiator(cdet, "Acrylic", CherenkovDetector::Upstream, isec,
0530                                         (G4LogicalVolume*)(0x2), 0, surface, filterThickness / mm);
0531           radiator->SetAlternativeMaterialName(filterMatName.c_str());
0532         }
0533       } //for isec
0534     }
0535 #endif
0536   }
0537 
0538 #ifdef WITH_IRT1_SUPPORT
0539   // radiator material names
0540   desc.add(Constant("DRICH_aerogel_material", aerogelMat.ptr()->GetName(), "string"));
0541   desc.add(Constant("DRICH_airgap_material", airgapMat.ptr()->GetName(), "string"));
0542   desc.add(Constant("DRICH_filter_material", filterMat.ptr()->GetName(), "string"));
0543   desc.add(Constant("DRICH_gasvol_material", gasvolMat.ptr()->GetName(), "string"));
0544 #endif
0545 #ifdef WITH_IRT2_SUPPORT
0546   // [0,0]: have neither access to G4VSolid nor to G4Material; IRT code does not care; fine;
0547   auto pd = new IRT2::CherenkovPhotonDetector(0, 0);
0548 
0549   // FIXME: '0' stands for the unknown (and irrelevant) G4LogicalVolume;
0550   geometry->AddPhotonDetector(cdet, 0, pd);
0551 
0552   // Cannot access GEANT shapes in the reconstruction code -> store this value;
0553   pd->SetActiveAreaSize(pssSide / mm);
0554 
0555   // FIXME: calculate it properly later; see S13361-3050NE-08 specs;
0556   pd->SetGeometricEfficiency(0.74);
0557 #endif
0558 
0559   // SECTOR LOOP //////////////////////////////////////////////////////////////////////
0560   for (int isec = 0; isec < nSectors; isec++) {
0561     // debugging filters, limiting the number of sectors
0562     if (debugSector && isec != 0)
0563       continue;
0564 
0565     // sector rotation about z axis
0566     RotationZ sectorRotation((isec + 0.5) * 2 * M_PI / nSectors);
0567     std::string secName = "sec" + std::to_string(isec);
0568 
0569     // BUILD MIRRORS ====================================================================
0570 
0571     // mirror positioning attributes
0572     // - sensor sphere center, w.r.t. IP
0573     double zS = sensorSphCenterZ + vesselZmin;
0574     double xS = sensorSphCenterX;
0575     // - distance between IP and mirror back plane
0576     double b = vesselZmax - mirrorBackplane;
0577     // - desired focal region: sensor sphere center, offset by focus-tune (z,x) parameters
0578     double zF = zS + focusTuneZ;
0579     double xF = xS + focusTuneX;
0580 
0581     // determine the mirror that focuses the IP to this desired region
0582     /* - uses point-to-point focusing to derive spherical mirror center
0583      *   `(mirrorCenterZ,mirrorCenterX)` and radius `mirrorRadius` for given
0584      *   image point coordinates `(zF,xF)` and `b`, defined as the z-distance
0585      *   between the object (IP) and the mirror surface
0586      * - all coordinates are specified w.r.t. the object point (IP)
0587      */
0588     double mirrorCenterZ = b * zF / (2 * b - zF);
0589     double mirrorCenterX = b * xF / (2 * b - zF);
0590     double mirrorRadius  = b - mirrorCenterZ;
0591 
0592     // translate mirror center to be w.r.t vessel front plane
0593     mirrorCenterZ -= vesselZmin;
0594 
0595     // spherical mirror patch cuts and rotation
0596     double mirrorThetaRot = std::asin(mirrorCenterX / mirrorRadius);
0597     double mirrorTheta1   = mirrorThetaRot - std::asin((mirrorCenterX - mirrorRmin) / mirrorRadius);
0598     double mirrorTheta2   = mirrorThetaRot + std::asin((mirrorRmax - mirrorCenterX) / mirrorRadius);
0599 
0600     // if debugging, draw full sphere
0601     if (debugMirror) {
0602       mirrorTheta1 = 0;
0603       mirrorTheta2 = M_PI;
0604     }
0605 
0606     // solid : create sphere at origin, with specified angular limits;
0607     // phi limits are increased to fill gaps (overlaps are cut away later)
0608     Sphere mirrorSolid1(mirrorRadius, mirrorRadius + mirrorThickness, mirrorTheta1, mirrorTheta2,
0609                         -40 * degree, 40 * degree);
0610 
0611     // mirror placement transformation (note: transformations are in reverse order)
0612     auto mirrorPos = Position(mirrorCenterX, 0., mirrorCenterZ) + originFront;
0613     auto mirrorPlacement(
0614         Translation3D(mirrorPos) * // re-center to specified position
0615         RotationY(-mirrorThetaRot) // rotate about vertical axis, to be within vessel radial walls
0616     );
0617 
0618     // cut overlaps with other sectors using "pie slice" wedges, to the extent specified
0619     // by `mirrorPhiw`
0620     Tube pieSlice(0.01 * cm, vesselRmax2, tankLength / 2.0, -mirrorPhiw / 2.0, mirrorPhiw / 2.0);
0621     IntersectionSolid mirrorSolid2(pieSlice, mirrorSolid1, mirrorPlacement);
0622 
0623     // mirror volume, attributes, and placement
0624     Volume mirrorVol(detName + "_mirror_" + secName, mirrorSolid2, mirrorMat);
0625     mirrorVol.setVisAttributes(mirrorVis);
0626     auto mirrorSectorPlacement = Transform3D(sectorRotation); // rotate about beam axis to sector
0627     auto mirrorPV              = gasvolVol.placeVolume(mirrorVol, mirrorSectorPlacement);
0628 
0629     // properties
0630     DetElement mirrorDE(det, "mirror_de_" + secName, isec);
0631     mirrorDE.setPlacement(mirrorPV);
0632     SkinSurface mirrorSkin(desc, mirrorDE, "mirror_optical_surface_" + secName, mirrorSurf,
0633                            mirrorVol);
0634     mirrorSkin.isValid();
0635 
0636 #if defined(WITH_IRT2_SUPPORT) || defined(WITH_IRT1_SUPPORT)
0637     // reconstruction constants (w.r.t. IP)
0638     // - access sector center after `sectorRotation`
0639     auto mirrorFinalPlacement = mirrorSectorPlacement * mirrorPlacement;
0640     auto mirrorFinalCenter    = vesselPos + mirrorFinalPlacement.Translation().Vect();
0641 #endif
0642 
0643 #ifdef WITH_IRT2_SUPPORT
0644     {
0645       // NB: default is concave, which is fine;
0646       msurface = new SphericalSurface(
0647           (1 / mm) * TVector3(mirrorFinalCenter.x(), mirrorFinalCenter.y(), mirrorFinalCenter.z()),
0648           (1 / mm) * mirrorRadius);
0649       mboundary = new IRT2::OpticalBoundary(cdet->GetRadiator("GasVolume"), msurface, false);
0650       // Need to store it in a separate call, see a comment in CherenkovDetector.h;
0651       cdet->StoreOpticalBoundary(mboundary);
0652 
0653       // Assign gas volume rear surface (this mirror) by hand;
0654       cdet->GetContainerVolume()->m_Borders[isec].second = msurface;
0655     }
0656 #endif
0657 #ifdef WITH_IRT1_SUPPORT
0658     desc.add(Constant("DRICH_mirror_center_x_" + secName, std::to_string(mirrorFinalCenter.x())));
0659     desc.add(Constant("DRICH_mirror_center_y_" + secName, std::to_string(mirrorFinalCenter.y())));
0660     desc.add(Constant("DRICH_mirror_center_z_" + secName, std::to_string(mirrorFinalCenter.z())));
0661     if (isec == 0)
0662       desc.add(Constant("DRICH_mirror_radius", std::to_string(mirrorRadius)));
0663 #endif
0664 
0665     // BUILD SENSORS ====================================================================
0666 
0667     // if debugging sphere properties, restrict number of sensors drawn
0668     if (debugSensors) {
0669       pssSide = 2 * M_PI * sensorSphRadius / 64;
0670     }
0671 
0672     // reconstruction constants
0673     auto sensorSphPos         = Position(sensorSphCenterX, 0., sensorSphCenterZ) + originFront;
0674     auto sensorSphFinalCenter = sectorRotation * Position(xS, 0.0, zS);
0675 #ifdef WITH_IRT1_SUPPORT
0676     desc.add(
0677         Constant("DRICH_sensor_sph_center_x_" + secName, std::to_string(sensorSphFinalCenter.x())));
0678     desc.add(
0679         Constant("DRICH_sensor_sph_center_y_" + secName, std::to_string(sensorSphFinalCenter.y())));
0680     desc.add(
0681         Constant("DRICH_sensor_sph_center_z_" + secName, std::to_string(sensorSphFinalCenter.z())));
0682     if (isec == 0)
0683       desc.add(Constant("DRICH_sensor_sph_radius", std::to_string(sensorSphRadius)));
0684 #endif
0685 
0686     // SENSOR MODULE LOOP ------------------------
0687     /* ALGORITHM: generate sphere of positions
0688      * - NOTE: there are two coordinate systems here:
0689      *   - "global" the main EPIC coordinate system
0690      *   - "generator" (vars end in `Gen`) is a local coordinate system for
0691      *     generating points on a sphere; it is related to the global system by
0692      *     a rotation; we do this so the "patch" (subset of generated
0693      *     positions) of sensors we choose to build is near the equator, where
0694      *     point distribution is more uniform
0695      * - PROCEDURE: loop over `thetaGen`, with subloop over `phiGen`, each divided evenly
0696      *   - the number of points to generate depends how many PDUs
0697      *     can fit within each ring of constant `thetaGen` or `phiGen`
0698      *   - we divide the relevant circumference by the PDU size, and this
0699      *     number is allowed to be a fraction, because likely we don't care about
0700      *     generating a full sphere and don't mind a "seam" at the overlap point
0701      *   - if we pick a patch of the sphere near the equator, and not near
0702      *     the poles or seam, the sensor distribution will appear uniform
0703      */
0704 
0705     // initialize PDU number for this sector
0706     int ipdu = 0;
0707 
0708     // calculate PDU pitch: the distance between two adjacent PDUs
0709     double pduPitch = pduNumSensors * resinSide + (pduNumSensors + 1) * pduSensorGap + pduGap;
0710 
0711     // BUILD PDU ASSEMBLY TEMPLATE -------------------------------------------------------
0712     /* All PDUs within a sector have identical internal geometry; only their placement in
0713      * the gas volume differs (each PDU points to a different spot on the sensor sphere).
0714      * By building one Assembly template per sector and placing it N times, TGeo can share
0715      * the same volume objects for all PDU copies, reducing the unique TGeoVolume count,
0716      * the associated memory footprint, and the cost of geometry lookups.
0717      *
0718      * PhysVolID allocation:
0719      *  - "sipm"   is set on the pssVol placement inside sensorAssembly (varies per slot)
0720      *  - "sector" and "pdu" are set on the pduAssembly placement in gasvolVol (varies per PDU)
0721      * The VolumeManager accumulates IDs from all levels when encoding a hit's cellID, so the
0722      * full {sector, pdu, sipm} triple is correctly reconstructed at simulation time.
0723      */
0724 
0725     /* begin building sensors and PDUs, where:
0726      * - sensor assembly: collection of all objects for a single SiPM
0727      * - photodetector unit (PDU) assembly: matrix of SiPMs with services
0728      *   - coordinate system: the "origin" of the assembly will be the center of the
0729      *     outermost surface of the photosensitive surface (pss)
0730      *     - reconstruction can access the sensor surface position from the sensor
0731      *       assembly origin, which will ultimately have coordinates w.r.t. to the IP after
0732      *       placement in the dRICH vessel
0733      *     - the pss is segmented into SiPM pixels; gaps between the pixels
0734      *       are accounted for in reconstruction, and each pixel reads out as a unique `cellID`
0735      *     - `cellID` to position conversion will give pixel centroids within the pss volume,
0736      *       (not exactly at the pss surface, but rather in the center of the pss volume,
0737      *       so keep in mind the very small offset)
0738      */
0739 
0740     // photosensitive surface (pss) and resin solids
0741     Box pssSolid(pssSide / 2., pssSide / 2., pssThickness / 2.);
0742     Box resinSolid(resinSide / 2., resinSide / 2., resinThickness / 2.);
0743 
0744     // embed pss solid in resin solid, by subtracting `pssSolid` from `resinSolid`
0745     SubtractionSolid resinSolidEmbedded(
0746         resinSolid, pssSolid,
0747         Transform3D(Translation3D(0., 0., (resinThickness - pssThickness) / 2.)));
0748 
0749     /* NOTE:
0750      * Here we could add gaps (size=`DRICH_pixel_gap`) between the pixels
0751      * as additional resin volumes, but this would require several more
0752      * iterative boolean operations, which may cause significant
0753      * performance slow downs in the simulation. Alternatively, one can
0754      * create a pixel gap mask with several disjoint, thin `Box` volumes
0755      * just outside the pss surface (no booleans required), but this
0756      * would amount to a very large number of additional volumes. Instead,
0757      * we have decided to apply pixel gap masking to the digitization
0758      * algorithm, downstream in reconstruction.
0759      */
0760 
0761     // pss and resin volumes (one per sector, shared across all PDU copies)
0762     Volume pssVol(detName + "_pss_" + secName, pssSolid, pssMat);
0763     Volume resinVol(detName + "_resin_" + secName, resinSolidEmbedded, resinMat);
0764     pssVol.setVisAttributes(pssVis);
0765     resinVol.setVisAttributes(resinVis);
0766 
0767     // sensitivity
0768     if (!debugOptics || debugOpticsMode == 3)
0769       pssVol.setSensitiveDetector(sens);
0770 
0771     // sensor optical surface: applied once to pssVol; covers all instances in all PDU copies
0772     if (!debugOptics || debugOpticsMode == 3) {
0773       SkinSurface pssSkin(desc, det, "sensor_optical_surface_" + secName, pssSurf, pssVol);
0774       pssSkin.isValid();
0775     }
0776 
0777     // PDU sensor grid parameters
0778     double pduSensorPitch     = resinSide + pduSensorGap;
0779     double pduSensorOffsetMax = pduSensorPitch * (pduNumSensors - 1) / 2.0;
0780 
0781     // placement transformations within sensorAssembly (identical for all sensors)
0782     // - set assembly origin to pss outermost surface centroid
0783     auto pssPlacement   = Transform3D(Translation3D(0., 0., -pssThickness / 2.0));
0784     auto resinPlacement = Transform3D(Translation3D(0., 0., -resinThickness / 2.0));
0785 
0786     // PDU assembly template for this sector
0787     Assembly pduAssembly(detName + "_pdu_" + secName);
0788 
0789     // fill sensor grid: one sensorAssembly per (sensorIx, sensorIy) slot
0790     // sensorPVs[isipm] holds the pssPV for that slot, used for DetElement creation later
0791     std::vector<PlacedVolume> sensorPVs;
0792     for (int sensorIx = 0; sensorIx < pduNumSensors; sensorIx++) {
0793       for (int sensorIy = 0; sensorIy < pduNumSensors; sensorIy++) {
0794         int isipm               = sensorIx * pduNumSensors + sensorIy;
0795         double pduSensorOffsetX = sensorIx * pduSensorPitch - pduSensorOffsetMax;
0796         double pduSensorOffsetY = sensorIy * pduSensorPitch - pduSensorOffsetMax;
0797         Assembly sensorAssembly(detName + "_sensor_" + secName + "_" + std::to_string(isipm));
0798         auto pssPV = sensorAssembly.placeVolume(pssVol, pssPlacement);
0799         sensorAssembly.placeVolume(resinVol, resinPlacement);
0800         pssPV.addPhysVolID("sipm", isipm); // "sector" and "pdu" go on the pduAssembly placement
0801         pduAssembly.placeVolume(
0802             sensorAssembly, Transform3D(Translation3D(pduSensorOffsetX, pduSensorOffsetY, 0.0)));
0803         sensorPVs.push_back(pssPV);
0804       }
0805     } // end PDU SiPM matrix template loop
0806 
0807     // front service volumes (one per sector, placed once in pduAssembly template)
0808     Transform3D frontServiceTransformation = Transform3D(Translation3D(0., 0., -resinThickness));
0809     for (xml::Collection_t serviceElem(pduElem.child(_Unicode(frontservices)), _Unicode(service));
0810          serviceElem; ++serviceElem) {
0811       auto serviceName      = serviceElem.attr<std::string>(_Unicode(name));
0812       auto serviceSide      = serviceElem.attr<double>(_Unicode(side));
0813       auto serviceThickness = serviceElem.attr<double>(_Unicode(thickness));
0814       auto serviceMat       = desc.material(serviceElem.attr<std::string>(_Unicode(material)));
0815       auto serviceVis       = desc.visAttributes(serviceElem.attr<std::string>(_Unicode(vis)));
0816       Box serviceSolid(serviceSide / 2.0, serviceSide / 2.0, serviceThickness / 2.0);
0817       Volume serviceVol(detName + "_" + serviceName + "_" + secName, serviceSolid, serviceMat);
0818       serviceVol.setVisAttributes(serviceVis);
0819       frontServiceTransformation =
0820           Transform3D(Translation3D(0., 0., -serviceThickness / 2.0)) * frontServiceTransformation;
0821       pduAssembly.placeVolume(serviceVol, frontServiceTransformation);
0822       frontServiceTransformation =
0823           Transform3D(Translation3D(0., 0., -serviceThickness / 2.0)) * frontServiceTransformation;
0824     }
0825 
0826     // circuit board volumes (one per sector)
0827     auto boardsElem = pduElem.child(_Unicode(boards));
0828     auto boardsMat  = desc.material(boardsElem.attr<std::string>(_Unicode(material)));
0829     auto boardsVis  = desc.visAttributes(boardsElem.attr<std::string>(_Unicode(vis)));
0830     Transform3D backServiceTransformation;
0831     for (xml::Collection_t boardElem(boardsElem, _Unicode(board)); boardElem; ++boardElem) {
0832       auto boardName      = boardElem.attr<std::string>(_Unicode(name));
0833       auto boardWidth     = boardElem.attr<double>(_Unicode(width));
0834       auto boardLength    = boardElem.attr<double>(_Unicode(length));
0835       auto boardThickness = boardElem.attr<double>(_Unicode(thickness));
0836       auto boardOffset    = boardElem.attr<double>(_Unicode(offset));
0837       Box boardSolid(boardWidth / 2.0, boardThickness / 2.0, boardLength / 2.0);
0838       Volume boardVol(detName + "_" + boardName + "+" + secName, boardSolid, boardsMat);
0839       boardVol.setVisAttributes(boardsVis);
0840       auto boardTransformation =
0841           Translation3D(0., boardOffset, -boardLength / 2.0) * frontServiceTransformation;
0842       pduAssembly.placeVolume(boardVol, boardTransformation);
0843       if (boardName == "RDO")
0844         backServiceTransformation =
0845             Translation3D(0., 0., -boardLength) * frontServiceTransformation;
0846     }
0847 
0848     // back service volumes (one per sector)
0849     for (xml::Collection_t serviceElem(pduElem.child(_Unicode(backservices)), _Unicode(service));
0850          serviceElem; ++serviceElem) {
0851       auto serviceName      = serviceElem.attr<std::string>(_Unicode(name));
0852       auto serviceSide      = serviceElem.attr<double>(_Unicode(side));
0853       auto serviceThickness = serviceElem.attr<double>(_Unicode(thickness));
0854       auto serviceMat       = desc.material(serviceElem.attr<std::string>(_Unicode(material)));
0855       auto serviceVis       = desc.visAttributes(serviceElem.attr<std::string>(_Unicode(vis)));
0856       Box serviceSolid(serviceSide / 2.0, serviceSide / 2.0, serviceThickness / 2.0);
0857       Volume serviceVol(detName + "_" + serviceName + "_" + secName, serviceSolid, serviceMat);
0858       serviceVol.setVisAttributes(serviceVis);
0859       backServiceTransformation =
0860           Transform3D(Translation3D(0., 0., -serviceThickness / 2.0)) * backServiceTransformation;
0861       pduAssembly.placeVolume(serviceVol, backServiceTransformation);
0862       backServiceTransformation =
0863           Transform3D(Translation3D(0., 0., -serviceThickness / 2.0)) * backServiceTransformation;
0864     }
0865 
0866     // PDU PLACEMENT LOOP ----------------------------------------------------------------
0867 
0868     // thetaGen loop: iterate less than "0.5 circumference / sensor size" times
0869     double nTheta = M_PI * sensorSphRadius / pduPitch;
0870     for (int t = 0; t < (int)(nTheta + 0.5); t++) {
0871       double thetaGen = t / ((double)nTheta) * M_PI;
0872 
0873       // phiGen loop: iterate less than "circumference at this latitude / sensor size" times
0874       double nPhi = 2 * M_PI * sensorSphRadius * std::sin(thetaGen) / pduPitch;
0875       for (int p = 0; p < (int)(nPhi + 0.5); p++) {
0876         double phiGen = p / ((double)nPhi) * 2 * M_PI - M_PI; // shift to [-pi,pi]
0877 
0878         // determine global phi and theta
0879         // - convert {radius,thetaGen,phiGen} -> {xGen,yGen,zGen}
0880         double xGen = sensorSphRadius * std::sin(thetaGen) * std::cos(phiGen);
0881         double yGen = sensorSphRadius * std::sin(thetaGen) * std::sin(phiGen);
0882         double zGen = sensorSphRadius * std::cos(thetaGen);
0883         // - convert {xGen,yGen,zGen} -> global {x,y,z} via rotation
0884         double x = zGen;
0885         double y = xGen;
0886         double z = yGen;
0887 
0888         // shift global coordinates so we can apply spherical patch cuts
0889         double zCheck   = z + sensorSphCenterZ;
0890         double xCheck   = x + sensorSphCenterX;
0891         double yCheck   = y;
0892         double rCheck   = std::hypot(xCheck, yCheck);
0893         double phiCheck = std::atan2(yCheck, xCheck);
0894 
0895         // patch cut
0896         bool patchCut = std::fabs(phiCheck) < sensorSphPatchPhiw && zCheck > sensorSphPatchZmin &&
0897                         rCheck > sensorSphPatchRmin && rCheck < sensorSphPatchRmax;
0898         if (debugSensors)
0899           patchCut = std::fabs(phiCheck) < sensorSphPatchPhiw;
0900         if (patchCut) {
0901 
0902           // PDU placement definition: describe how to place a PDU on the sphere
0903           /* - transformations operate on global coordinates; the corresponding
0904            *   generator coordinates are provided in the comments
0905            * - transformations are applied in reverse order
0906            * - the `pduAssembly` origin is at the active surface; in other words, this origin
0907            *   should be placed on the sensor sphere surface by `pduAssemblyPlacement`
0908            */
0909           // clang-format off
0910           auto pduAssemblyPlacement =
0911             sectorRotation *                         // rotate about beam axis to sector
0912             Translation3D(sensorSphPos) *            // move sphere to reference position
0913             RotationX(phiGen) *                      // rotate about `zGen`
0914             RotationZ(thetaGen) *                    // rotate about `yGen`
0915             Translation3D(sensorSphRadius, 0., 0.) * // push radially to spherical surface
0916             RotationY(M_PI / 2) *                    // rotate sensor to be compatible with generator coords
0917             RotationZ(-M_PI / 2);                    // correction for readout segmentation mapping
0918           // clang-format on
0919 
0920           // place shared pduAssembly template; set "sector" and "pdu" physVolIDs on this placement
0921           auto pduPV = gasvolVol.placeVolume(pduAssembly, pduAssemblyPlacement);
0922           pduPV.addPhysVolID("sector", isec).addPhysVolID("pdu", ipdu);
0923 
0924           // per-sensor DetElement creation and parameter storage
0925           // (these depend on pduAssemblyPlacement which varies per PDU position)
0926           int isipm = 0;
0927           for (int sensorIx = 0; sensorIx < pduNumSensors; sensorIx++) {
0928             for (int sensorIy = 0; sensorIy < pduNumSensors; sensorIy++) {
0929 
0930               auto pduSensorOffsetX = sensorIx * pduSensorPitch - pduSensorOffsetMax;
0931               auto pduSensorOffsetY = sensorIy * pduSensorPitch - pduSensorOffsetMax;
0932               auto sensorAssemblyPlacement =
0933                   Transform3D(Translation3D(pduSensorOffsetX, pduSensorOffsetY, 0.0));
0934 
0935               // sensor readout: build sensorID from known (sector, pdu, sipm) values
0936               // NOTE: physVolIDs are split across levels: {sector,pdu} on pduPV, {sipm} on sensorPVs[isipm]
0937               auto sensorID = encodeSensorID(std::vector<std::pair<std::string, int>>{
0938                   {"sector", isec}, {"pdu", ipdu}, {"sipm", isipm}});
0939               //printf("@S@ %d vs %lu\n", isec, (sensorID >> 8) & 0x7);
0940               std::string sensorIDname =
0941                   secName + "_pdu" + std::to_string(ipdu) + "_sipm" + std::to_string(isipm);
0942               DetElement pssDE(det, "sensor_de_" + sensorIDname, sensorID);
0943               // Use pduPV (the concrete placed PDU instance) rather than the shared template
0944               // sensor node: pduPV is unique per PDU placement and carries the {sector,pdu}
0945               // physVolIDs, making DE-based lookups unambiguous. The individual sensor
0946               // position within the PDU is available via the VariantParameters below.
0947               pssDE.setPlacement(pduPV);
0948 
0949               // obtain some parameters useful for optics, so we don't have to figure them out downstream
0950               // - sensor position: the centroid of the active SURFACE of the `pss`
0951               auto pduOrigin = ROOT::Math::XYZPoint(0, 0, 0);
0952               auto sensorPos = Translation3D(vesselPos) * // position of vessel in world
0953                                pduAssemblyPlacement *     // position of PDU in vessel
0954                                sensorAssemblyPlacement *  // position of SiPM in PDU
0955                                pduOrigin;
0956               auto pduPos    = Translation3D(vesselPos) * // position of vessel in world
0957                                pduAssemblyPlacement *     // position of PDU in vessel
0958                                pduOrigin;
0959               // - sensor surface basis: the orientation of the sensor surface
0960               //   NOTE: all sensors of a single PDU have the same surface orientation, but to avoid
0961               //         loss of generality downstream, define the basis for each sensor
0962               auto normVector = [pduAssemblyPlacement](Direction n) {
0963                 return pduAssemblyPlacement * n;
0964               };
0965               auto sensorNormX = normVector(Direction{
0966                   -1., //1.,
0967                   0.,
0968                   0.,
0969               });
0970               auto sensorNormY = normVector(Direction{
0971                   0.,
0972                   -1., //1.,
0973                   0.,
0974               });
0975 
0976               // geometry tests
0977               /* - to help ensure the optics geometry is correctly interpreted by the reconstruction,
0978                *   we do a few checks here
0979                * - if any changes break these tests, the determination of `sensorPos`,
0980                *   `sensorNormX`, `sensorNormY` is wrong and/or the tests need to be updated
0981                */
0982               // - test: check if the sensor position is on the sensor sphere (corrected for PDU matrix offset)
0983               auto distActual   = std::sqrt((sensorPos - sensorSphFinalCenter).Mag2());
0984               auto distExpected = std::hypot(pduSensorOffsetX, pduSensorOffsetY, sensorSphRadius);
0985               auto testOnSphere = distActual - distExpected;
0986               if (std::abs(testOnSphere) > 1e-6) {
0987                 printout(ERROR, "DRICH_geo", "sensor %s failed on-sphere test; testOnSphere=%f",
0988                          sensorIDname.c_str(), testOnSphere);
0989                 throw std::runtime_error("dRICH sensor position test failed");
0990               }
0991               // - test: check the orientation
0992               Direction radialDir =
0993                   Direction(pduPos) - sensorSphFinalCenter;      // sensor sphere radius direction
0994               auto sensorNormZ = sensorNormX.Cross(sensorNormY); // sensor surface normal
0995               auto testOrtho =
0996                   sensorNormX.Dot(sensorNormY); // zero, if x and y vectors are orthogonal
0997               auto testRadial =
0998                   radialDir.Cross(sensorNormZ)
0999                       .Mag2(); // zero, if surface normal is parallel to radial direction
1000               auto testDirection = radialDir.Dot(
1001                   sensorNormZ); // positive, if radial direction == sensor normal direction (outward)
1002               if (std::abs(testOrtho) > 1e-6 || std::abs(testRadial) > 1e-6 || testDirection <= 0) {
1003                 printout(ERROR, "DRICH_geo", "sensor %s failed orientation test",
1004                          sensorIDname.c_str());
1005                 printout(ERROR, "DRICH_geo", "  testOrtho     = %f; should be zero", testOrtho);
1006                 printout(ERROR, "DRICH_geo", "  testRadial    = %f; should be zero", testRadial);
1007                 printout(ERROR, "DRICH_geo", "  testDirection = %f; should be positive",
1008                          testDirection);
1009                 throw std::runtime_error("dRICH sensor orientation test failed");
1010               }
1011 
1012 #ifdef WITH_IRT2_SUPPORT
1013               {
1014                 // SiPM panel surface;
1015                 auto surface = new FlatSurface(
1016                     (1 / mm) * TVector3(sensorPos.x(), sensorPos.y(), sensorPos.z()),
1017                     TVector3(sensorNormX.x(), sensorNormX.y(), sensorNormX.z()),
1018                     TVector3(sensorNormY.x(), sensorNormY.y(), sensorNormY.z()));
1019 
1020                 // Wipe out sector bits; FIXME: it seems this is not really needed?;
1021                 auto irt = pd->AllocateIRT(isec, sensorID & sector_mask);
1022 
1023                 // Aerogel and acrylic;
1024                 if (cdet->m_OpticalBoundaries[CherenkovDetector::Upstream].find(isec) !=
1025                     cdet->m_OpticalBoundaries[CherenkovDetector::Upstream].end())
1026                   for (auto boundary : cdet->m_OpticalBoundaries[CherenkovDetector::Upstream][isec])
1027                     irt->AddOpticalBoundary(boundary);
1028 
1029                 // Mirror;
1030                 irt->AddOpticalBoundary(mboundary);
1031 
1032                 // FIXME: eventually there should be a quartz window defined as part of the
1033                 // cdet->m_OpticalBoundaries[CherenkovDetector::Downstream] boundaries;
1034 
1035                 // Terminate the optical path;
1036                 pd->AddItselfToOpticalBoundaries(irt, surface);
1037               }
1038 #endif
1039 #ifdef WITH_IRT1_SUPPORT
1040               // add these optics parameters to this sensor's parameter map
1041               auto pssVarMap = pssDE.extension<VariantParameters>(false);
1042               if (pssVarMap == nullptr) {
1043                 pssVarMap = new VariantParameters();
1044                 pssDE.addExtension<VariantParameters>(pssVarMap);
1045               }
1046 
1047               auto addVecToMap = [pssVarMap](std::string key, auto vec) {
1048                 pssVarMap->set<double>(key + "_x", vec.x());
1049                 pssVarMap->set<double>(key + "_y", vec.y());
1050                 pssVarMap->set<double>(key + "_z", vec.z());
1051               };
1052               addVecToMap("pos", sensorPos);
1053               addVecToMap("normX", sensorNormX);
1054               addVecToMap("normY", sensorNormY);
1055               printout(DEBUG, "DRICH_geo", "sensor %s:", sensorIDname.c_str());
1056               for (auto kv : pssVarMap->variantParameters)
1057                 printout(DEBUG, "DRICH_geo", "    %s: %f", kv.first.c_str(),
1058                          pssVarMap->get<double>(kv.first));
1059 #endif
1060               printout(DEBUG, "DRICH_geo", "sensor %s:", sensorIDname.c_str());
1061 
1062               // increment SIPM number
1063               isipm++;
1064             }
1065           } // end per-PDU sensor DetElement loop
1066 
1067           // increment PDU number
1068           ipdu++;
1069 
1070         } // end patch cuts
1071       } // end phiGen loop
1072     } // end thetaGen loop
1073 
1074     // END SENSOR MODULE LOOP ------------------------
1075 
1076 #ifdef WITH_IRT1_SUPPORT
1077     // add constant for access to the number of PDUs per sector
1078     if (isec == 0)
1079       desc.add(Constant("DRICH_num_pdus", std::to_string(ipdu)));
1080     else if (ipdu != desc.constant<int>("DRICH_num_pdus"))
1081       printout(WARNING, "DRICH_geo", "number of PDUs is not the same for each sector");
1082 #endif
1083 
1084   } // END SECTOR LOOP //////////////////////////
1085 
1086   return det;
1087 }
1088 
1089 // clang-format off
1090 DECLARE_DETELEMENT(epic_DRICH, createDetector)
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!