epic/src/DRICH_geo.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022, 2023 Christopher Dilks, Junhuai Xu
0003
0004 //==========================================================================
0005 //  dRICH: Dual Ring Imaging Cherenkov Detector
0006 //--------------------------------------------------------------------------
0007 //
0008 // Author: Christopher Dilks (Duke University)
0009 //
0010 // - Design Adapted from Standalone Fun4all and GEMC implementations
0011 //   [ Evaristo Cisbani, Cristiano Fanelli, Alessio Del Dotto, et al. ]
0012 //
0013 //==========================================================================
0014
0015 #include "DD4hep/DetFactoryHelper.h"
0016 #include "DD4hep/OpticalSurfaces.h"
0017 #include "DD4hep/Printout.h"
0018 #include "DDRec/DetectorData.h"
0019 #include "DDRec/Surface.h"
0020
0021 #include <XML/Helper.h>
0022
0023 using namespace dd4hep;
0024 using namespace dd4hep::rec;
0025
0026 // create the detector
0027 static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetector sens) {
0028
0029   xml::DetElement detElem       = handle;
0030   std::string detName           = detElem.nameStr();
0031   int detID                     = detElem.id();
0032   xml::Component dims           = detElem.dimensions();
0033   OpticalSurfaceManager surfMgr = desc.surfaceManager();
0034   DetElement det(detName, detID);
0035   sens.setType("tracker");
0036
0037   // attributes, from compact file =============================================
0038   // - vessel
0039   auto vesselZmin      = dims.attr<double>(_Unicode(zmin));
0040   auto vesselLength    = dims.attr<double>(_Unicode(length));
0041   auto vesselRmin0     = dims.attr<double>(_Unicode(rmin0));
0042   auto vesselRmin1     = dims.attr<double>(_Unicode(rmin1));
0043   auto vesselRmax0     = dims.attr<double>(_Unicode(rmax0));
0044   auto vesselRmax1     = dims.attr<double>(_Unicode(rmax1));
0045   auto vesselRmax2     = dims.attr<double>(_Unicode(rmax2));
0046   auto snoutLength     = dims.attr<double>(_Unicode(snout_length));
0047   auto nSectors        = dims.attr<int>(_Unicode(nsectors));
0048   auto wallThickness   = dims.attr<double>(_Unicode(wall_thickness));
0049   auto windowThickness = dims.attr<double>(_Unicode(window_thickness));
0050   auto vesselMat       = desc.material(detElem.attr<std::string>(_Unicode(material)));
0051   auto gasvolMat       = desc.material(detElem.attr<std::string>(_Unicode(gas)));
0052   auto vesselVis       = desc.visAttributes(detElem.attr<std::string>(_Unicode(vis_vessel)));
0053   auto gasvolVis       = desc.visAttributes(detElem.attr<std::string>(_Unicode(vis_gas)));
0054   // - radiator (applies to aerogel and filter)
0055   auto radiatorElem       = detElem.child(_Unicode(radiator));
0056   auto radiatorRmin       = radiatorElem.attr<double>(_Unicode(rmin));
0057   auto radiatorRmax       = radiatorElem.attr<double>(_Unicode(rmax));
0058   auto radiatorPitch      = radiatorElem.attr<double>(_Unicode(pitch));
0059   auto radiatorFrontplane = radiatorElem.attr<double>(_Unicode(frontplane));
0060   // - aerogel
0061   auto aerogelElem      = radiatorElem.child(_Unicode(aerogel));
0062   auto aerogelMat       = desc.material(aerogelElem.attr<std::string>(_Unicode(material)));
0063   auto aerogelVis       = desc.visAttributes(aerogelElem.attr<std::string>(_Unicode(vis)));
0064   auto aerogelThickness = aerogelElem.attr<double>(_Unicode(thickness));
0065   // - filter
0066   auto filterElem      = radiatorElem.child(_Unicode(filter));
0067   auto filterMat       = desc.material(filterElem.attr<std::string>(_Unicode(material)));
0068   auto filterVis       = desc.visAttributes(filterElem.attr<std::string>(_Unicode(vis)));
0069   auto filterThickness = filterElem.attr<double>(_Unicode(thickness));
0070   // - airgap between filter and aerogel
0071   auto airgapElem      = radiatorElem.child(_Unicode(airgap));
0072   auto airgapMat       = desc.material(airgapElem.attr<std::string>(_Unicode(material)));
0073   auto airgapVis       = desc.visAttributes(airgapElem.attr<std::string>(_Unicode(vis)));
0074   auto airgapThickness = airgapElem.attr<double>(_Unicode(thickness));
0075   // - mirror
0076   auto mirrorElem      = detElem.child(_Unicode(mirror));
0077   auto mirrorMat       = desc.material(mirrorElem.attr<std::string>(_Unicode(material)));
0078   auto mirrorVis       = desc.visAttributes(mirrorElem.attr<std::string>(_Unicode(vis)));
0079   auto mirrorSurf      = surfMgr.opticalSurface(mirrorElem.attr<std::string>(_Unicode(surface)));
0080   auto mirrorBackplane = mirrorElem.attr<double>(_Unicode(backplane));
0081   auto mirrorThickness = mirrorElem.attr<double>(_Unicode(thickness));
0082   auto mirrorRmin      = mirrorElem.attr<double>(_Unicode(rmin));
0083   auto mirrorRmax      = mirrorElem.attr<double>(_Unicode(rmax));
0084   auto mirrorPhiw      = mirrorElem.attr<double>(_Unicode(phiw));
0085   auto focusTuneZ      = mirrorElem.attr<double>(_Unicode(focus_tune_z));
0086   auto focusTuneX      = mirrorElem.attr<double>(_Unicode(focus_tune_x));
0087   // - sensorboxes
0088   auto sensorboxLength = desc.constant<double>("DRICH_sensorbox_length");
0089   auto sensorboxRmin   = desc.constant<double>("DRICH_sensorbox_rmin");
0090   auto sensorboxRmax   = desc.constant<double>("DRICH_sensorbox_rmax");
0091   auto sensorboxDphi   = desc.constant<double>("DRICH_sensorbox_dphi");
0092   // - sensor photosensitive surface (pss)
0093   auto pssElem      = detElem.child(_Unicode(sensors)).child(_Unicode(pss));
0094   auto pssMat       = desc.material(pssElem.attr<std::string>(_Unicode(material)));
0095   auto pssVis       = desc.visAttributes(pssElem.attr<std::string>(_Unicode(vis)));
0096   auto pssSurf      = surfMgr.opticalSurface(pssElem.attr<std::string>(_Unicode(surface)));
0097   auto pssSide      = pssElem.attr<double>(_Unicode(side));
0098   auto pssThickness = pssElem.attr<double>(_Unicode(thickness));
0099   // - sensor resin
0100   auto resinElem      = detElem.child(_Unicode(sensors)).child(_Unicode(resin));
0101   auto resinMat       = desc.material(resinElem.attr<std::string>(_Unicode(material)));
0102   auto resinVis       = desc.visAttributes(resinElem.attr<std::string>(_Unicode(vis)));
0103   auto resinSide      = resinElem.attr<double>(_Unicode(side));
0104   auto resinThickness = resinElem.attr<double>(_Unicode(thickness));
0105   // - photodetector unit (PDU)
0106   auto pduElem       = detElem.child(_Unicode(sensors)).child(_Unicode(pdu));
0107   auto pduNumSensors = desc.constant<int>("DRICH_pdu_num_sensors");
0108   auto pduSensorGap  = desc.constant<double>("DRICH_pdu_sensor_gap");
0109   auto pduGap        = desc.constant<double>("DRICH_pdu_gap");
0110   // - sensor sphere
0111   auto sensorSphElem    = detElem.child(_Unicode(sensors)).child(_Unicode(sphere));
0112   auto sensorSphRadius  = sensorSphElem.attr<double>(_Unicode(radius));
0113   auto sensorSphCenterX = sensorSphElem.attr<double>(_Unicode(centerx));
0114   auto sensorSphCenterZ = sensorSphElem.attr<double>(_Unicode(centerz));
0115   // - sensor sphere patch cuts
0116   auto sensorSphPatchElem = detElem.child(_Unicode(sensors)).child(_Unicode(sphericalpatch));
0117   auto sensorSphPatchPhiw = sensorSphPatchElem.attr<double>(_Unicode(phiw));
0118   auto sensorSphPatchRmin = sensorSphPatchElem.attr<double>(_Unicode(rmin));
0119   auto sensorSphPatchRmax = sensorSphPatchElem.attr<double>(_Unicode(rmax));
0120   auto sensorSphPatchZmin = sensorSphPatchElem.attr<double>(_Unicode(zmin));
0121   // - sensor readout
0122   auto readoutName = detElem.attr<std::string>(_Unicode(readout));
0123   // - settings and switches
0124   auto debugOpticsMode = desc.constant<int>("DRICH_debug_optics");
0125   bool debugSector     = desc.constant<int>("DRICH_debug_sector") == 1;
0126   bool debugMirror     = desc.constant<int>("DRICH_debug_mirror") == 1;
0127   bool debugSensors    = desc.constant<int>("DRICH_debug_sensors") == 1;
0128
0129   // if debugging optics, override some settings
0130   bool debugOptics = debugOpticsMode > 0;
0131   if (debugOptics) {
0132     printout(WARNING, "DRICH_geo", "DEBUGGING DRICH OPTICS");
0133     switch (debugOpticsMode) {
0134     case 1:
0135       vesselMat = aerogelMat = filterMat = pssMat = gasvolMat = desc.material("VacuumOptical");
0136       break;
0137     case 2:
0138       vesselMat = aerogelMat = filterMat = pssMat = desc.material("VacuumOptical");
0139       break;
0140     case 3:
0141       vesselMat = aerogelMat = filterMat = gasvolMat = desc.material("VacuumOptical");
0142       break;
0143     default:
0144       printout(FATAL, "DRICH_geo", "UNKNOWN debugOpticsMode");
0145       return det;
0146     }
0147   }
0148
0149   // if debugging anything, draw only one sector and adjust visibility
0150   if (debugOptics || debugMirror || debugSensors)
0151     debugSector = true;
0152   if (debugSector)
0153     gasvolVis = vesselVis = desc.invisible();
0154
0155   // readout coder <-> unique sensor ID
0156   /* - `sensorIDfields` is a list of readout fields used to specify a unique sensor ID
0157    * - `cellMask` is defined such that a hit's `cellID & cellMask` is the corresponding sensor's unique ID
0158    */
0159   std::vector<std::string> sensorIDfields = {"pdu", "sipm", "sector"};
0160   const auto& readoutCoder                = *desc.readout(readoutName).idSpec().decoder();
0161   // determine `cellMask` based on `sensorIDfields`
0162   uint64_t cellMask = 0;
0163   for (const auto& idField : sensorIDfields)
0164     cellMask |= readoutCoder[idField].mask();
0165   desc.add(Constant("DRICH_cell_mask", std::to_string(cellMask)));
0166   // create a unique sensor ID from a sensor's PlacedVolume::volIDs
0167   auto encodeSensorID = [&readoutCoder](auto ids) {
0168     uint64_t enc = 0;
0169     for (const auto& [idField, idValue] : ids)
0170       enc |= uint64_t(idValue) << readoutCoder[idField].offset();
0171     return enc;
0172   };
0173
0174   // BUILD VESSEL ====================================================================
0175   /* - `vessel`: aluminum enclosure, the mother volume of the dRICH
0176    * - `gasvol`: gas volume, which fills `vessel`; all other volumes defined below
0177    *   are children of `gasvol`
0178    * - the dRICH vessel geometry has two regions: the snout refers to the conic region
0179    *   in the front, housing the aerogel, while the tank refers to the cylindrical
0180    *   region, housing the rest of the detector components
0181    */
0182
0183   // derived attributes
0184   double tankLength = vesselLength - snoutLength;
0185   double vesselZmax = vesselZmin + vesselLength;
0186
0187   // snout solids
0188   double boreDelta  = vesselRmin1 - vesselRmin0;
0189   double snoutDelta = vesselRmax1 - vesselRmax0;
0190   Cone vesselSnout(snoutLength / 2.0, vesselRmin0, vesselRmax0,
0191                    vesselRmin0 + boreDelta * snoutLength / vesselLength, vesselRmax1);
0192   Cone gasvolSnout(
0193       /* note: `gasvolSnout` extends a bit into the tank, so it touches `gasvolTank`
0194        * - the extension distance is equal to the tank `windowThickness`, so the
0195        *   length of `gasvolSnout` == length of `vesselSnout`
0196        * - the extension backplane radius is calculated using similar triangles
0197        */
0198       snoutLength / 2.0, vesselRmin0 + wallThickness, vesselRmax0 - wallThickness,
0199       vesselRmin0 + boreDelta * (snoutLength - windowThickness) / vesselLength + wallThickness,
0200       vesselRmax1 - wallThickness + windowThickness * (vesselRmax1 - vesselRmax0) / snoutLength);
0201
0202   // tank solids:
0203   // - inner: cone along beamline
0204   // - outer: cone to back of sensor box, then fixed radius cylinder
0205   Polycone vesselTank(
0206       0, 2 * M_PI,
0207       /* rmin */
0208       {vesselSnout.rMin2(),
0209        std::lerp(vesselSnout.rMin2(), vesselRmin1, (sensorboxLength - snoutLength) / tankLength),
0210        vesselRmin1},
0211       /* rmax */ {vesselSnout.rMax2(), vesselRmax2, vesselRmax2},
0212       /* z    */
0213       {-tankLength / 2.0, -tankLength / 2.0 + sensorboxLength - snoutLength, tankLength / 2.0});
0214   Polycone gasvolTank(
0215       0, 2 * M_PI,
0216       /* rmin */
0217       {gasvolSnout.rMin2(),
0218        std::lerp(gasvolSnout.rMin2(), vesselRmin1 + wallThickness,
0219                  (sensorboxLength - snoutLength) / tankLength),
0220        vesselRmin1 + wallThickness},
0221       /* rmax */ {gasvolSnout.rMax2(), vesselRmax2 - wallThickness, vesselRmax2 - wallThickness},
0222       /* z    */
0223       {-tankLength / 2.0 + windowThickness,
0224        -tankLength / 2.0 + windowThickness + sensorboxLength - snoutLength,
0225        tankLength / 2.0 - windowThickness});
0226
0227   // sensorbox solids
0228   double dphi = atan2(wallThickness, sensorboxRmax); // thickness only correct at Rmax
0229   Tube vesselSensorboxTube(sensorboxRmin, sensorboxRmax, sensorboxLength / 2., -sensorboxDphi / 2.,
0230                            sensorboxDphi / 2.);
0231   Tube gasvolSensorboxTube(sensorboxRmin + wallThickness, sensorboxRmax - wallThickness,
0232                            sensorboxLength / 2., -sensorboxDphi / 2. + dphi,
0233                            sensorboxDphi / 2. - dphi);
0234
0235   // union: snout + tank
0236   UnionSolid vesselUnion(vesselTank, vesselSnout, Position(0., 0., -vesselLength / 2.));
0237   UnionSolid gasvolUnion(gasvolTank, gasvolSnout,
0238                          Position(0., 0., -vesselLength / 2. + windowThickness));
0239
0240   // union: add sensorboxes for all sectors
0241   for (int isec = 0; isec < nSectors; isec++) {
0242     RotationZ sectorRotation((isec + 0.5) * 2 * M_PI / nSectors);
0243     vesselUnion = UnionSolid(
0244         vesselUnion, vesselSensorboxTube,
0245         Transform3D(sectorRotation, Position(0., 0., -(snoutLength + sensorboxLength - 0.6) / 2.)));
0246     gasvolUnion = UnionSolid(
0247         gasvolUnion, gasvolSensorboxTube,
0248         Transform3D(sectorRotation,
0249                     Position(0., 0., -(snoutLength + sensorboxLength) / 2. + windowThickness)));
0250   }
0251
0252   //  extra solids for `debugOptics` only
0253   Box vesselBox(1001, 1001, 1001);
0254   Box gasvolBox(1000, 1000, 1000);
0255
0256   // choose vessel and gasvol solids (depending on `debugOpticsMode` (0=disabled))
0257   Solid vesselSolid, gasvolSolid;
0258   switch (debugOpticsMode) {
0259   case 0:
0260     vesselSolid = vesselUnion;
0261     gasvolSolid = gasvolUnion;
0262     break; // `!debugOptics`
0263   case 1:
0264   case 3:
0265     vesselSolid = vesselBox;
0266     gasvolSolid = gasvolBox;
0267     break;
0268   case 2:
0269     vesselSolid = vesselBox;
0270     gasvolSolid = gasvolUnion;
0271     break;
0272   }
0273
0274   // volumes
0275   Volume vesselVol(detName, vesselSolid, vesselMat);
0276   Volume gasvolVol(detName + "_gas", gasvolSolid, gasvolMat);
0277   vesselVol.setVisAttributes(vesselVis);
0278   gasvolVol.setVisAttributes(gasvolVis);
0279
0280   // reference positions
0281   // - the vessel is created such that the center of the cylindrical tank volume
0282   //   coincides with the origin; this is called the "origin position" of the vessel
0283   // - when the vessel (and its children volumes) is placed, it is translated in
0284   //   the z-direction to be in the proper EPIC-integration location
0285   // - these reference positions are for the frontplane and backplane of the vessel,
0286   //   with respect to the vessel origin position
0287   auto originFront = Position(0., 0., -tankLength / 2.0 - snoutLength);
0288   // auto originBack  = Position(0., 0., tankLength / 2.0);
0289   auto vesselPos = Position(0, 0, vesselZmin) - originFront;
0290
0291   // place gas volume
0292   PlacedVolume gasvolPV = vesselVol.placeVolume(gasvolVol, Position(0, 0, 0));
0293   DetElement gasvolDE(det, "gasvol_de", 0);
0294   gasvolDE.setPlacement(gasvolPV);
0295
0296   // place mother volume (vessel)
0297   Volume motherVol      = desc.pickMotherVolume(det);
0298   PlacedVolume vesselPV = motherVol.placeVolume(vesselVol, vesselPos);
0299   vesselPV.addPhysVolID("system", detID);
0300   det.setPlacement(vesselPV);
0301
0302   // BUILD RADIATOR ====================================================================
0303
0304   // solid and volume: create aerogel and filter
0305   Cone aerogelSolid(aerogelThickness / 2, radiatorRmin, radiatorRmax,
0306                     radiatorRmin + boreDelta * aerogelThickness / vesselLength,
0307                     radiatorRmax + snoutDelta * aerogelThickness / snoutLength);
0308   Cone airgapSolid(airgapThickness / 2, radiatorRmin + boreDelta * aerogelThickness / vesselLength,
0309                    radiatorRmax + snoutDelta * aerogelThickness / snoutLength,
0310                    radiatorRmin + boreDelta * (aerogelThickness + airgapThickness) / vesselLength,
0311                    radiatorRmax + snoutDelta * (aerogelThickness + airgapThickness) / snoutLength);
0312   Cone filterSolid(
0313       filterThickness / 2,
0314       radiatorRmin + boreDelta * (aerogelThickness + airgapThickness) / vesselLength,
0315       radiatorRmax + snoutDelta * (aerogelThickness + airgapThickness) / snoutLength,
0316       radiatorRmin +
0317           boreDelta * (aerogelThickness + airgapThickness + filterThickness) / vesselLength,
0318       radiatorRmax +
0319           snoutDelta * (aerogelThickness + airgapThickness + filterThickness) / snoutLength);
0320
0321   Volume aerogelVol(detName + "_aerogel", aerogelSolid, aerogelMat);
0322   Volume airgapVol(detName + "_airgap", airgapSolid, airgapMat);
0323   Volume filterVol(detName + "_filter", filterSolid, filterMat);
0324   aerogelVol.setVisAttributes(aerogelVis);
0325   airgapVol.setVisAttributes(airgapVis);
0326   filterVol.setVisAttributes(filterVis);
0327
0328   // aerogel placement and surface properties
0329   // TODO [low-priority]: define skin properties for aerogel and filter
0330   // FIXME: radiatorPitch might not be working correctly (not yet used)
0331   auto radiatorPos = Position(0., 0., radiatorFrontplane + 0.5 * aerogelThickness) + originFront;
0332   auto aerogelPlacement = Translation3D(radiatorPos) * // re-center to originFront
0333                           RotationY(radiatorPitch);    // change polar angle to specified pitch
0334   auto aerogelPV = gasvolVol.placeVolume(aerogelVol, aerogelPlacement);
0335   DetElement aerogelDE(det, "aerogel_de", 0);
0336   aerogelDE.setPlacement(aerogelPV);
0337   // SkinSurface aerogelSkin(desc, aerogelDE, "mirror_optical_surface", aerogelSurf, aerogelVol);
0338   // aerogelSkin.isValid();
0339
0340   // airgap and filter placement and surface properties
0341   if (!debugOptics) {
0342
0343     auto airgapPlacement =
0344         Translation3D(radiatorPos) * // re-center to originFront
0345         RotationY(radiatorPitch) *   // change polar angle
0346         Translation3D(0., 0.,
0347                       (aerogelThickness + airgapThickness) / 2.); // move to aerogel backplane
0348     auto airgapPV = gasvolVol.placeVolume(airgapVol, airgapPlacement);
0349     DetElement airgapDE(det, "airgap_de", 0);
0350     airgapDE.setPlacement(airgapPV);
0351
0352     auto filterPlacement =
0353         Translation3D(0., 0., airgapThickness) * // add an air gap
0354         Translation3D(radiatorPos) *             // re-center to originFront
0355         RotationY(radiatorPitch) *               // change polar angle
0356         Translation3D(0., 0.,
0357                       (aerogelThickness + filterThickness) / 2.); // move to aerogel backplane
0358     auto filterPV = gasvolVol.placeVolume(filterVol, filterPlacement);
0359     DetElement filterDE(det, "filter_de", 0);
0360     filterDE.setPlacement(filterPV);
0361     // SkinSurface filterSkin(desc, filterDE, "mirror_optical_surface", filterSurf, filterVol);
0362     // filterSkin.isValid();
0363
0364     // radiator z-positions (w.r.t. IP); only needed downstream if !debugOptics
0365     double aerogelZpos = vesselPos.z() + aerogelPV.position().z();
0366     double airgapZpos  = vesselPos.z() + airgapPV.position().z();
0367     double filterZpos  = vesselPos.z() + filterPV.position().z();
0368     desc.add(Constant("DRICH_aerogel_zpos", std::to_string(aerogelZpos)));
0369     desc.add(Constant("DRICH_airgap_zpos", std::to_string(airgapZpos)));
0370     desc.add(Constant("DRICH_filter_zpos", std::to_string(filterZpos)));
0371   }
0372
0373   // radiator material names
0374   desc.add(Constant("DRICH_aerogel_material", aerogelMat.ptr()->GetName(), "string"));
0375   desc.add(Constant("DRICH_airgap_material", airgapMat.ptr()->GetName(), "string"));
0376   desc.add(Constant("DRICH_filter_material", filterMat.ptr()->GetName(), "string"));
0377   desc.add(Constant("DRICH_gasvol_material", gasvolMat.ptr()->GetName(), "string"));
0378
0379   // SECTOR LOOP //////////////////////////////////////////////////////////////////////
0380   for (int isec = 0; isec < nSectors; isec++) {
0381
0382     // debugging filters, limiting the number of sectors
0383     if (debugSector && isec != 0)
0384       continue;
0385
0386     // sector rotation about z axis
0387     RotationZ sectorRotation((isec + 0.5) * 2 * M_PI / nSectors);
0388     std::string secName = "sec" + std::to_string(isec);
0389
0390     // BUILD MIRRORS ====================================================================
0391
0392     // mirror positioning attributes
0393     // - sensor sphere center, w.r.t. IP
0394     double zS = sensorSphCenterZ + vesselZmin;
0395     double xS = sensorSphCenterX;
0396     // - distance between IP and mirror back plane
0397     double b = vesselZmax - mirrorBackplane;
0398     // - desired focal region: sensor sphere center, offset by focus-tune (z,x) parameters
0399     double zF = zS + focusTuneZ;
0400     double xF = xS + focusTuneX;
0401
0402     // determine the mirror that focuses the IP to this desired region
0403     /* - uses point-to-point focusing to derive spherical mirror center
0404      *   `(mirrorCenterZ,mirrorCenterX)` and radius `mirrorRadius` for given
0405      *   image point coordinates `(zF,xF)` and `b`, defined as the z-distance
0406      *   between the object (IP) and the mirror surface
0407      * - all coordinates are specified w.r.t. the object point (IP)
0408      */
0409     double mirrorCenterZ = b * zF / (2 * b - zF);
0410     double mirrorCenterX = b * xF / (2 * b - zF);
0411     double mirrorRadius  = b - mirrorCenterZ;
0412
0413     // translate mirror center to be w.r.t vessel front plane
0414     mirrorCenterZ -= vesselZmin;
0415
0416     // spherical mirror patch cuts and rotation
0417     double mirrorThetaRot = std::asin(mirrorCenterX / mirrorRadius);
0418     double mirrorTheta1   = mirrorThetaRot - std::asin((mirrorCenterX - mirrorRmin) / mirrorRadius);
0419     double mirrorTheta2   = mirrorThetaRot + std::asin((mirrorRmax - mirrorCenterX) / mirrorRadius);
0420
0421     // if debugging, draw full sphere
0422     if (debugMirror) {
0423       mirrorTheta1 = 0;
0424       mirrorTheta2 = M_PI; /*mirrorPhiw=2*M_PI;*/
0425     }
0426
0427     // solid : create sphere at origin, with specified angular limits;
0428     // phi limits are increased to fill gaps (overlaps are cut away later)
0429     Sphere mirrorSolid1(mirrorRadius, mirrorRadius + mirrorThickness, mirrorTheta1, mirrorTheta2,
0430                         -40 * degree, 40 * degree);
0431
0432     // mirror placement transformation (note: transformations are in reverse order)
0433     auto mirrorPos = Position(mirrorCenterX, 0., mirrorCenterZ) + originFront;
0434     auto mirrorPlacement(
0435         Translation3D(mirrorPos) * // re-center to specified position
0436         RotationY(-mirrorThetaRot) // rotate about vertical axis, to be within vessel radial walls
0437     );
0438
0439     // cut overlaps with other sectors using "pie slice" wedges, to the extent specified
0440     // by `mirrorPhiw`
0441     Tube pieSlice(0.01 * cm, vesselRmax2, tankLength / 2.0, -mirrorPhiw / 2.0, mirrorPhiw / 2.0);
0442     IntersectionSolid mirrorSolid2(pieSlice, mirrorSolid1, mirrorPlacement);
0443
0444     // mirror volume, attributes, and placement
0445     Volume mirrorVol(detName + "_mirror_" + secName, mirrorSolid2, mirrorMat);
0446     mirrorVol.setVisAttributes(mirrorVis);
0447     auto mirrorSectorPlacement = Transform3D(sectorRotation); // rotate about beam axis to sector
0448     auto mirrorPV              = gasvolVol.placeVolume(mirrorVol, mirrorSectorPlacement);
0449
0450     // properties
0451     DetElement mirrorDE(det, "mirror_de_" + secName, isec);
0452     mirrorDE.setPlacement(mirrorPV);
0453     SkinSurface mirrorSkin(desc, mirrorDE, "mirror_optical_surface_" + secName, mirrorSurf,
0454                            mirrorVol);
0455     mirrorSkin.isValid();
0456
0457     // reconstruction constants (w.r.t. IP)
0458     // - access sector center after `sectorRotation`
0459     auto mirrorFinalPlacement = mirrorSectorPlacement * mirrorPlacement;
0460     auto mirrorFinalCenter    = vesselPos + mirrorFinalPlacement.Translation().Vect();
0461     desc.add(Constant("DRICH_mirror_center_x_" + secName, std::to_string(mirrorFinalCenter.x())));
0462     desc.add(Constant("DRICH_mirror_center_y_" + secName, std::to_string(mirrorFinalCenter.y())));
0463     desc.add(Constant("DRICH_mirror_center_z_" + secName, std::to_string(mirrorFinalCenter.z())));
0464     if (isec == 0)
0465       desc.add(Constant("DRICH_mirror_radius", std::to_string(mirrorRadius)));
0466
0467     // BUILD SENSORS ====================================================================
0468
0469     // if debugging sphere properties, restrict number of sensors drawn
0470     if (debugSensors) {
0471       pssSide = 2 * M_PI * sensorSphRadius / 64;
0472     }
0473
0474     // reconstruction constants
0475     auto sensorSphPos         = Position(sensorSphCenterX, 0., sensorSphCenterZ) + originFront;
0476     auto sensorSphFinalCenter = sectorRotation * Position(xS, 0.0, zS);
0477     desc.add(
0478         Constant("DRICH_sensor_sph_center_x_" + secName, std::to_string(sensorSphFinalCenter.x())));
0479     desc.add(
0480         Constant("DRICH_sensor_sph_center_y_" + secName, std::to_string(sensorSphFinalCenter.y())));
0481     desc.add(
0482         Constant("DRICH_sensor_sph_center_z_" + secName, std::to_string(sensorSphFinalCenter.z())));
0483     if (isec == 0)
0484       desc.add(Constant("DRICH_sensor_sph_radius", std::to_string(sensorSphRadius)));
0485
0486     // SENSOR MODULE LOOP ------------------------
0487     /* ALGORITHM: generate sphere of positions
0488      * - NOTE: there are two coordinate systems here:
0489      *   - "global" the main EPIC coordinate system
0490      *   - "generator" (vars end in `Gen`) is a local coordinate system for
0491      *     generating points on a sphere; it is related to the global system by
0492      *     a rotation; we do this so the "patch" (subset of generated
0493      *     positions) of sensors we choose to build is near the equator, where
0494      *     point distribution is more uniform
0495      * - PROCEDURE: loop over `thetaGen`, with subloop over `phiGen`, each divided evenly
0496      *   - the number of points to generate depends how many PDUs
0497      *     can fit within each ring of constant `thetaGen` or `phiGen`
0498      *   - we divide the relevant circumference by the PDU size, and this
0499      *     number is allowed to be a fraction, because likely we don't care about
0500      *     generating a full sphere and don't mind a "seam" at the overlap point
0501      *   - if we pick a patch of the sphere near the equator, and not near
0502      *     the poles or seam, the sensor distribution will appear uniform
0503      */
0504
0505     // initialize PDU number for this sector
0506     int ipdu = 0;
0507
0508     // calculate PDU pitch: the distance between two adjacent PDUs
0509     double pduPitch = pduNumSensors * resinSide + (pduNumSensors + 1) * pduSensorGap + pduGap;
0510
0511     // thetaGen loop: iterate less than "0.5 circumference / sensor size" times
0512     double nTheta = M_PI * sensorSphRadius / pduPitch;
0513     for (int t = 0; t < (int)(nTheta + 0.5); t++) {
0514       double thetaGen = t / ((double)nTheta) * M_PI;
0515
0516       // phiGen loop: iterate less than "circumference at this latitude / sensor size" times
0517       double nPhi = 2 * M_PI * sensorSphRadius * std::sin(thetaGen) / pduPitch;
0518       for (int p = 0; p < (int)(nPhi + 0.5); p++) {
0519         double phiGen = p / ((double)nPhi) * 2 * M_PI - M_PI; // shift to [-pi,pi]
0520
0521         // determine global phi and theta
0522         // - convert {radius,thetaGen,phiGen} -> {xGen,yGen,zGen}
0523         double xGen = sensorSphRadius * std::sin(thetaGen) * std::cos(phiGen);
0524         double yGen = sensorSphRadius * std::sin(thetaGen) * std::sin(phiGen);
0525         double zGen = sensorSphRadius * std::cos(thetaGen);
0526         // - convert {xGen,yGen,zGen} -> global {x,y,z} via rotation
0527         double x = zGen;
0528         double y = xGen;
0529         double z = yGen;
0530         // - convert global {x,y,z} -> global {phi,theta}
0531         // double phi   = std::atan2(y, x);
0532         // double theta = std::acos(z / sensorSphRadius);
0533
0534         // shift global coordinates so we can apply spherical patch cuts
0535         double zCheck   = z + sensorSphCenterZ;
0536         double xCheck   = x + sensorSphCenterX;
0537         double yCheck   = y;
0538         double rCheck   = std::hypot(xCheck, yCheck);
0539         double phiCheck = std::atan2(yCheck, xCheck);
0540
0541         // patch cut
0542         bool patchCut = std::fabs(phiCheck) < sensorSphPatchPhiw && zCheck > sensorSphPatchZmin &&
0543                         rCheck > sensorSphPatchRmin && rCheck < sensorSphPatchRmax;
0544         if (debugSensors)
0545           patchCut = std::fabs(phiCheck) < sensorSphPatchPhiw;
0546         if (patchCut) {
0547
0548           /* begin building sensors and PDUs, where:
0549            * - sensor assembly: collection of all objects for a single SiPM
0550            * - photodetector unit (PDU) assembly: matrix of SiPMs with services
0551            *   - coordinate system: the "origin" of the assembly will be the center of the
0552            *     outermost surface of the photosensitive surface (pss)
0553            *     - reconstruction can access the sensor surface position from the sensor
0554            *       assembly origin, which will ultimately have coordinates w.r.t. to the IP after
0555            *       placement in the dRICH vessel
0556            *     - the pss is segmented into SiPM pixels; gaps between the pixels
0557            *       are accounted for in reconstruction, and each pixel reads out as a unique `cellID`
0558            *     - `cellID` to postion conversion will give pixel centroids within the pss volume,
0559            *       (not exactly at the pss surface, but rather in the center of the pss volume,
0560            *       so keep in mind the very small offset)
0561            */
0562
0563           // photosensitive surface (pss) and resin solids
0564           Box pssSolid(pssSide / 2., pssSide / 2., pssThickness / 2.);
0565           Box resinSolid(resinSide / 2., resinSide / 2., resinThickness / 2.);
0566
0567           // embed pss solid in resin solid, by subtracting `pssSolid` from `resinSolid`
0568           SubtractionSolid resinSolidEmbedded(
0569               resinSolid, pssSolid,
0570               Transform3D(Translation3D(0., 0., (resinThickness - pssThickness) / 2.)));
0571
0572           /* NOTE:
0573            * Here we could add gaps (size=`DRICH_pixel_gap`) between the pixels
0574            * as additional resin volumes, but this would require several more
0575            * iterative boolean operations, which may cause significant
0576            * performance slow downs in the simulation. Alternatively, one can
0577            * create a pixel gap mask with several disjoint, thin `Box` volumes
0578            * just outside the pss surface (no booleans required), but this
0579            * would amount to a very large number of additional volumes. Instead,
0580            * we have decided to apply pixel gap masking to the digitization
0581            * algorithm, downstream in reconstruction.
0582            */
0583
0584           // pss and resin volumes
0585           Volume pssVol(detName + "_pss_" + secName, pssSolid, pssMat);
0586           Volume resinVol(detName + "_resin_" + secName, resinSolidEmbedded, resinMat);
0587           pssVol.setVisAttributes(pssVis);
0588           resinVol.setVisAttributes(resinVis);
0589
0590           // sensitivity
0591           if (!debugOptics || debugOpticsMode == 3)
0592             pssVol.setSensitiveDetector(sens);
0593
0594           // PDU placement definition: describe how to place a PDU on the sphere
0595           /* - transformations operate on global coordinates; the corresponding
0596            *   generator coordinates are provided in the comments
0597            * - transformations are applied in reverse order
0598            * - the `pduAssembly` origin is at the active surface; in other words, this origin
0599            *   should be placed on the sensor sphere surface by `pduAssemblyPlacement`
0600            */
0601           // clang-format off
0602           auto pduAssemblyPlacement =
0603             sectorRotation *                         // rotate about beam axis to sector
0604             Translation3D(sensorSphPos) *            // move sphere to reference position
0605             RotationX(phiGen) *                      // rotate about `zGen`
0606             RotationZ(thetaGen) *                    // rotate about `yGen`
0607             Translation3D(sensorSphRadius, 0., 0.) * // push radially to spherical surface
0608             RotationY(M_PI / 2) *                    // rotate sensor to be compatible with generator coords
0609             RotationZ(-M_PI / 2);                    // correction for readout segmentation mapping
0610           // clang-format on
0611
0612           // generate matrix of sensors and place them in `pduAssembly`
0613           Assembly pduAssembly(detName + "_pdu_" + secName);
0614           double pduSensorPitch     = resinSide + pduSensorGap;
0615           double pduSensorOffsetMax = pduSensorPitch * (pduNumSensors - 1) / 2.0;
0616           int isipm                 = 0;
0617           for (int sensorIx = 0; sensorIx < pduNumSensors; sensorIx++) {
0618             for (int sensorIy = 0; sensorIy < pduNumSensors; sensorIy++) {
0619
0620               Assembly sensorAssembly(detName + "_sensor_" + secName);
0621
0622               // placement transformations
0623               // - placement of objects in `sensorAssembly`
0624               auto pssPlacement   = Transform3D(Translation3D(
0625                   0., 0.,
0626                   -pssThickness / 2.0)); // set assembly origin to pss outermost surface centroid
0627               auto resinPlacement = Transform3D(Translation3D(0., 0., -resinThickness / 2.0));
0628               // - placement of a `sensorAssembly` in `pduAssembly`
0629               auto pduSensorOffsetX = sensorIx * pduSensorPitch - pduSensorOffsetMax;
0630               auto pduSensorOffsetY = sensorIy * pduSensorPitch - pduSensorOffsetMax;
0631               auto sensorAssemblyPlacement =
0632                   Transform3D(Translation3D(pduSensorOffsetX, pduSensorOffsetY, 0.0));
0633
0634               // placements
0635               auto pssPV = sensorAssembly.placeVolume(pssVol, pssPlacement);
0636               sensorAssembly.placeVolume(resinVol, resinPlacement);
0637               pduAssembly.placeVolume(sensorAssembly, sensorAssemblyPlacement);
0638
0639               // sensor readout // NOTE: follow `sensorIDfields`
0640               pssPV.addPhysVolID("sector", isec)
0641                   .addPhysVolID("pdu", ipdu)
0642                   .addPhysVolID("sipm", isipm);
0643
0644               // sensor DetElement
0645               auto sensorID = encodeSensorID(pssPV.volIDs());
0646               std::string sensorIDname =
0647                   secName + "_pdu" + std::to_string(ipdu) + "_sipm" + std::to_string(isipm);
0648               DetElement pssDE(det, "sensor_de_" + sensorIDname, sensorID);
0649               pssDE.setPlacement(pssPV);
0650
0651               // sensor surface properties
0652               if (!debugOptics || debugOpticsMode == 3) {
0653                 SkinSurface pssSkin(desc, pssDE, "sensor_optical_surface_" + sensorIDname, pssSurf,
0654                                     pssVol);
0655                 pssSkin.isValid();
0656               }
0657
0658               // obtain some parameters useful for optics, so we don't have to figure them out downstream
0659               // - sensor position: the centroid of the active SURFACE of the `pss`
0660               auto pduOrigin = ROOT::Math::XYZPoint(0, 0, 0);
0661               auto sensorPos = Translation3D(vesselPos) * // position of vessel in world
0662                                pduAssemblyPlacement *     // position of PDU in vessel
0663                                sensorAssemblyPlacement *  // position of SiPM in PDU
0664                                pduOrigin;
0665               auto pduPos = Translation3D(vesselPos) * // position of vessel in world
0666                             pduAssemblyPlacement *     // position of PDU in vessel
0667                             pduOrigin;
0668               // - sensor surface basis: the orientation of the sensor surface
0669               //   NOTE: all sensors of a single PDU have the same surface orientation, but to avoid
0670               //         loss of generality downstream, define the basis for each sensor
0671               auto normVector = [pduAssemblyPlacement](Direction n) {
0672                 return pduAssemblyPlacement * n;
0673               };
0674               auto sensorNormX = normVector(Direction{
0675                   1.,
0676                   0.,
0677                   0.,
0678               });
0679               auto sensorNormY = normVector(Direction{
0680                   0.,
0681                   1.,
0682                   0.,
0683               });
0684
0685               // geometry tests
0686               /* - to help ensure the optics geometry is correctly interpreted by the reconstruction,
0687                *   we do a few checks here
0688                * - if any changes break these tests, the determination of `sensorPos`,
0689                *   `sensorNormX`, `sensorNormY` is wrong and/or the tests need to be updated
0690                */
0691               // - test: check if the sensor position is on the sensor sphere (corrected for PDU matrix offset)
0692               auto distActual   = std::sqrt((sensorPos - sensorSphFinalCenter).Mag2());
0693               auto distExpected = std::hypot(pduSensorOffsetX, pduSensorOffsetY, sensorSphRadius);
0694               auto testOnSphere = distActual - distExpected;
0695               if (std::abs(testOnSphere) > 1e-6) {
0696                 printout(ERROR, "DRICH_geo", "sensor %s failed on-sphere test; testOnSphere=%f",
0697                          sensorIDname.c_str(), testOnSphere);
0698                 throw std::runtime_error("dRICH sensor position test failed");
0699               }
0700               // - test: check the orientation
0701               Direction radialDir =
0702                   Direction(pduPos) - sensorSphFinalCenter;      // sensor sphere radius direction
0703               auto sensorNormZ = sensorNormX.Cross(sensorNormY); // sensor surface normal
0704               auto testOrtho =
0705                   sensorNormX.Dot(sensorNormY); // zero, if x and y vectors are orthogonal
0706               auto testRadial =
0707                   radialDir.Cross(sensorNormZ)
0708                       .Mag2(); // zero, if surface normal is parallel to radial direction
0709               auto testDirection = radialDir.Dot(
0710                   sensorNormZ); // positive, if radial direction == sensor normal direction (outward)
0711               if (std::abs(testOrtho) > 1e-6 || std::abs(testRadial) > 1e-6 || testDirection <= 0) {
0712                 printout(ERROR, "DRICH_geo", "sensor %s failed orientation test",
0713                          sensorIDname.c_str());
0714                 printout(ERROR, "DRICH_geo", "  testOrtho     = %f; should be zero", testOrtho);
0715                 printout(ERROR, "DRICH_geo", "  testRadial    = %f; should be zero", testRadial);
0716                 printout(ERROR, "DRICH_geo", "  testDirection = %f; should be positive",
0717                          testDirection);
0718                 throw std::runtime_error("dRICH sensor orientation test failed");
0719               }
0720
0721               // add these optics parameters to this sensor's parameter map
0722               auto pssVarMap = pssDE.extension<VariantParameters>(false);
0723               if (pssVarMap == nullptr) {
0724                 pssVarMap = new VariantParameters();
0725                 pssDE.addExtension<VariantParameters>(pssVarMap);
0726               }
0727               auto addVecToMap = [pssVarMap](std::string key, auto vec) {
0728                 pssVarMap->set<double>(key + "_x", vec.x());
0729                 pssVarMap->set<double>(key + "_y", vec.y());
0730                 pssVarMap->set<double>(key + "_z", vec.z());
0731               };
0732               addVecToMap("pos", sensorPos);
0733               addVecToMap("normX", sensorNormX);
0734               addVecToMap("normY", sensorNormY);
0735               printout(DEBUG, "DRICH_geo", "sensor %s:", sensorIDname.c_str());
0736               for (auto kv : pssVarMap->variantParameters)
0737                 printout(DEBUG, "DRICH_geo", "    %s: %f", kv.first.c_str(),
0738                          pssVarMap->get<double>(kv.first));
0739
0740               // increment SIPM number
0741               isipm++;
0742             }
0743           } // end PDU SiPM matrix loop
0744
0745           // front service volumes
0746           Transform3D frontServiceTransformation =
0747               Transform3D(Translation3D(0., 0., -resinThickness));
0748           for (xml::Collection_t serviceElem(pduElem.child(_Unicode(frontservices)),
0749                                              _Unicode(service));
0750                serviceElem; ++serviceElem) {
0751             auto serviceName      = serviceElem.attr<std::string>(_Unicode(name));
0752             auto serviceSide      = serviceElem.attr<double>(_Unicode(side));
0753             auto serviceThickness = serviceElem.attr<double>(_Unicode(thickness));
0754             auto serviceMat = desc.material(serviceElem.attr<std::string>(_Unicode(material)));
0755             auto serviceVis = desc.visAttributes(serviceElem.attr<std::string>(_Unicode(vis)));
0756             Box serviceSolid(serviceSide / 2.0, serviceSide / 2.0, serviceThickness / 2.0);
0757             Volume serviceVol(detName + "_" + serviceName + "_" + secName, serviceSolid,
0758                               serviceMat);
0759             serviceVol.setVisAttributes(serviceVis);
0760             frontServiceTransformation =
0761                 Transform3D(Translation3D(0., 0., -serviceThickness / 2.0)) *
0762                 frontServiceTransformation;
0763             pduAssembly.placeVolume(serviceVol, frontServiceTransformation);
0764             frontServiceTransformation =
0765                 Transform3D(Translation3D(0., 0., -serviceThickness / 2.0)) *
0766                 frontServiceTransformation;
0767           }
0768
0769           // circuit board volumes
0770           auto boardsElem = pduElem.child(_Unicode(boards));
0771           auto boardsMat  = desc.material(boardsElem.attr<std::string>(_Unicode(material)));
0772           auto boardsVis  = desc.visAttributes(boardsElem.attr<std::string>(_Unicode(vis)));
0773           Transform3D backServiceTransformation;
0774           for (xml::Collection_t boardElem(boardsElem, _Unicode(board)); boardElem; ++boardElem) {
0775             auto boardName      = boardElem.attr<std::string>(_Unicode(name));
0776             auto boardWidth     = boardElem.attr<double>(_Unicode(width));
0777             auto boardLength    = boardElem.attr<double>(_Unicode(length));
0778             auto boardThickness = boardElem.attr<double>(_Unicode(thickness));
0779             auto boardOffset    = boardElem.attr<double>(_Unicode(offset));
0780             Box boardSolid(boardWidth / 2.0, boardThickness / 2.0, boardLength / 2.0);
0781             Volume boardVol(detName + "_" + boardName + "+" + secName, boardSolid, boardsMat);
0782             boardVol.setVisAttributes(boardsVis);
0783             auto boardTransformation =
0784                 Translation3D(0., boardOffset, -boardLength / 2.0) * frontServiceTransformation;
0785             pduAssembly.placeVolume(boardVol, boardTransformation);
0786             if (boardName == "RDO")
0787               backServiceTransformation =
0788                   Translation3D(0., 0., -boardLength) * frontServiceTransformation;
0789           }
0790
0791           // back service volumes
0792           for (xml::Collection_t serviceElem(pduElem.child(_Unicode(backservices)),
0793                                              _Unicode(service));
0794                serviceElem; ++serviceElem) {
0795             auto serviceName      = serviceElem.attr<std::string>(_Unicode(name));
0796             auto serviceSide      = serviceElem.attr<double>(_Unicode(side));
0797             auto serviceThickness = serviceElem.attr<double>(_Unicode(thickness));
0798             auto serviceMat = desc.material(serviceElem.attr<std::string>(_Unicode(material)));
0799             auto serviceVis = desc.visAttributes(serviceElem.attr<std::string>(_Unicode(vis)));
0800             Box serviceSolid(serviceSide / 2.0, serviceSide / 2.0, serviceThickness / 2.0);
0801             Volume serviceVol(detName + "_" + serviceName + "_" + secName, serviceSolid,
0802                               serviceMat);
0803             serviceVol.setVisAttributes(serviceVis);
0804             backServiceTransformation =
0805                 Transform3D(Translation3D(0., 0., -serviceThickness / 2.0)) *
0806                 backServiceTransformation;
0807             pduAssembly.placeVolume(serviceVol, backServiceTransformation);
0808             backServiceTransformation =
0809                 Transform3D(Translation3D(0., 0., -serviceThickness / 2.0)) *
0810                 backServiceTransformation;
0811           }
0812
0813           // place PDU assembly
0814           gasvolVol.placeVolume(pduAssembly, pduAssemblyPlacement);
0815
0816           // increment PDU number
0817           ipdu++;
0818
0819         } // end patch cuts
0820       } // end phiGen loop
0821     } // end thetaGen loop
0822
0823     // END SENSOR MODULE LOOP ------------------------
0824
0825     // add constant for access to the number of PDUs per sector
0826     if (isec == 0)
0827       desc.add(Constant("DRICH_num_pdus", std::to_string(ipdu)));
0828     else if (ipdu != desc.constant<int>("DRICH_num_pdus"))
0829       printout(WARNING, "DRICH_geo", "number of PDUs is not the same for each sector");
0830
0831   } // END SECTOR LOOP //////////////////////////
0832
0833   return det;
0834 }
0835
0836 // clang-format off
0837 DECLARE_DETELEMENT(epic_DRICH, createDetector)