epic/src/PFRICH_geo.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2023  Wenliang (Bill) Li, Alexander Kiselev, Karthik Suresh
0003
0004 //----------------------------------
0005 // pfRICH: Proximity Focusing RICH
0006 // Author: Wenliang (Bill) Li
0007 //
0008 // - Design Adapted from standalone Geant4 description by
0009 //   Alexander Kiselev and Chandradoy Chatterjee
0010 //----------------------------------
0011
0012 #include "DD4hep/DetFactoryHelper.h"
0013
0014 #include "TVector2.h"
0015
0016 using namespace dd4hep;
0017
0018 #ifdef WITH_IRT2_SUPPORT
0019 #include "IRT2/CherenkovDetectorCollection.h"
0020 #include "IRT2/ConicalSurface.h"
0021
0022 using namespace IRT2;
0023 #endif
0024
0025 // -------------------------------------------------------------------------------------
0026
0027 static UnionSolid FlangeCut(Detector& description, double length, double clearance) {
0028   // FIXME: not the most efficient way to recalculate them every time new,
0029   // but code is more readable, and overhead is negligible anyway;
0030   auto _FLANGE_EPIPE_DIAMETER_ = description.constant<double>("FLANGE_EPIPE_DIAMETER");
0031   auto _FLANGE_HPIPE_DIAMETER_ = description.constant<double>("FLANGE_HPIPE_DIAMETER");
0032   auto _FLANGE_HPIPE_OFFSET_   = description.constant<double>("FLANGE_HPIPE_OFFSET");
0033
0034   // A wedge bridging two cylinders;
0035   Tube eflange(0.0, _FLANGE_EPIPE_DIAMETER_ / 2 + clearance, length / 2);
0036   Tube hflange(0.0, _FLANGE_HPIPE_DIAMETER_ / 2 + clearance, length / 2);
0037
0038   double r0 = _FLANGE_EPIPE_DIAMETER_ / 2 + clearance;
0039   double r1 = _FLANGE_HPIPE_DIAMETER_ / 2 + clearance;
0040   double L  = _FLANGE_HPIPE_OFFSET_;
0041   double a  = r0 * L / (r0 - r1);
0042   double b  = r0 * r0 / a;
0043   double c  = r1 * (a - b) / r0;
0044
0045   // GEANT variables to define G4Trap;
0046   double pDz = length / 2, pTheta = 0.0, pPhi = 0.0, pDy1 = (a - b - c) / 2, pDy2 = pDy1;
0047   double pDx1 = sqrt(r0 * r0 - b * b), pDx2 = pDx1 * r1 / r0, pDx3 = pDx1, pDx4 = pDx2, pAlp1 = 0.0,
0048          pAlp2 = 0.0;
0049
0050   Trap wedge(pDz, pTheta, pPhi, pDy1, pDx1, pDx2, pAlp1, pDy2, pDx3, pDx4, pAlp2);
0051
0052   UnionSolid flange_shape(eflange, hflange, Position(-_FLANGE_HPIPE_OFFSET_, 0.0, 0.0));
0053   Rotation3D rZ(RotationZYX(M_PI / 2.0, 0.0, 0.0));
0054   Transform3D transform_flange(rZ, Position(-b - pDy1, 0.0, 0.0));
0055
0056   return UnionSolid(flange_shape, wedge, transform_flange);
0057 } // FlangeCut()
0058
0059 // -------------------------------------------------------------------------------------
0060
0061 static Ref_t createDetector(Detector& description, xml_h e, SensitiveDetector sens) {
0062   xml::DetElement detElem = e;
0063   int det_id              = detElem.id();
0064
0065   std::string det_name = detElem.nameStr();
0066   Material air         = description.air();
0067
0068   DetElement sdet(det_name, det_id);
0069
0070   sens.setType("tracker");
0071   description.invisible();
0072
0073   std::string detName = detElem.nameStr();
0074
0075   int id = detElem.hasAttr(_U(id)) ? detElem.id() : 0;
0076
0077   // Start optical configuration if needed;
0078 #ifdef WITH_IRT2_SUPPORT
0079   auto geometry = CherenkovDetectorCollection::Instance();
0080   auto cdet     = geometry->AddNewDetector(detName.c_str());
0081 #endif
0082
0083   OpticalSurfaceManager surfMgr = description.surfaceManager();
0084
0085   auto aerogelElem      = detElem.child(_Unicode(aerogel));
0086   auto aerogelThickness = aerogelElem.attr<double>(_Unicode(thickness));
0087
0088   auto filterElem      = detElem.child(_Unicode(filter));
0089   auto filterThickness = filterElem.attr<double>(_Unicode(thickness));
0090
0091   // readout coder <-> unique sensor ID
0092   /* - `sensorIDfields` is a list of readout fields used to specify a unique sensor ID
0093      * - `cellMask` is defined such that a hit's `cellID & cellMask` is the corresponding sensor's unique ID
0094      * - this redundant generalization is for future flexibility, and consistency with dRICH
0095      */
0096   std::vector<std::string> sensorIDfields = {"hrppd"};
0097   auto readoutName                        = detElem.attr<std::string>(_Unicode(readout));
0098   const auto& readoutCoder                = *description.readout(readoutName).idSpec().decoder();
0099   // determine `cellMask` based on `sensorIDfields`
0100   uint64_t cellMask = 0;
0101   for (const auto& idField : sensorIDfields)
0102     cellMask |= readoutCoder[idField].mask();
0103   description.add(Constant("PFRICH_cell_mask", std::to_string(cellMask)));
0104 #ifdef WITH_IRT2_SUPPORT
0105   // Do not mind to store it twice;
0106   cdet->SetReadoutCellMask(cellMask);
0107 #endif
0108   // create a unique sensor ID from a sensor's PlacedVolume::volIDs
0109   auto encodeSensorID = [&readoutCoder](auto ids) {
0110     uint64_t enc = 0;
0111     for (const auto& [idField, idValue] : ids)
0112       enc |= uint64_t(idValue) << readoutCoder[idField].offset();
0113     return enc;
0114   };
0115
0116 #ifdef WITH_IRT2_SUPPORT
0117   const double sign = -1.0;
0118 #endif
0119
0120   //
0121   // Fiducial volume; will be split into three subvolumes along the beam line: front wall, gas volume
0122   // and sensor compartment; FIXME: called 'air' but presently filled with nitrogen;
0123   //
0124   auto _FIDUCIAL_VOLUME_LENGTH_ = description.constant<double>("FIDUCIAL_VOLUME_LENGTH");
0125   auto _VESSEL_OUTER_RADIUS_    = description.constant<double>("VESSEL_OUTER_RADIUS");
0126
0127   Tube pfRICH_air_volume(0.0, _VESSEL_OUTER_RADIUS_, _FIDUCIAL_VOLUME_LENGTH_ / 2);
0128
0129   auto _FLANGE_CLEARANCE_ = description.constant<double>("FLANGE_CLEARANCE");
0130   SubtractionSolid pfRICH_volume_shape(
0131       pfRICH_air_volume,
0132       FlangeCut(description, _FIDUCIAL_VOLUME_LENGTH_ + 1 * mm, _FLANGE_CLEARANCE_));
0133   auto vesselGasName = detElem.attr<std::string>(_Unicode(gas));
0134   auto vesselGas     = description.material(vesselGasName);
0135   Volume pfRICH_volume(detName, pfRICH_volume_shape, vesselGas);
0136   auto gasvolVis = description.visAttributes(detElem.attr<std::string>(_Unicode(vis_gas)));
0137   pfRICH_volume.setVisAttributes(gasvolVis);
0138   pfRICH_volume.setRegion(description, detElem.regionStr());
0139
0140   Volume mother                 = description.pickMotherVolume(sdet);
0141   auto _FIDUCIAL_VOLUME_OFFSET_ = description.constant<double>("FIDUCIAL_VOLUME_OFFSET");
0142   Transform3D transform(RotationZYX(0, M_PI, 0), Position(0, 0, _FIDUCIAL_VOLUME_OFFSET_));
0143   PlacedVolume pv = mother.placeVolume(pfRICH_volume, transform);
0144   if (id != 0)
0145     pv.addPhysVolID("system", id);
0146   sdet.setPlacement(pv);
0147
0148   //
0149   // Gas volume
0150   //
0151   auto _VESSEL_FRONT_SIDE_THICKNESS_ = description.constant<double>("VESSEL_FRONT_SIDE_THICKNESS");
0152   auto _SENSOR_AREA_LENGTH_          = description.constant<double>("SENSOR_AREA_LENGTH");
0153   auto _VESSEL_OUTER_WALL_THICKNESS_ = description.constant<double>("VESSEL_OUTER_WALL_THICKNESS");
0154
0155   // FIXME: do it better later;
0156 #ifdef WITH_IRT2_SUPPORT
0157   double fvOffset = fabs(_FIDUCIAL_VOLUME_OFFSET_);
0158 #endif
0159
0160   double gas_volume_length =
0161       _FIDUCIAL_VOLUME_LENGTH_ - _VESSEL_FRONT_SIDE_THICKNESS_ - _SENSOR_AREA_LENGTH_;
0162   double gas_volume_radius = _VESSEL_OUTER_RADIUS_ - _VESSEL_OUTER_WALL_THICKNESS_;
0163   double gas_volume_offset = -(_SENSOR_AREA_LENGTH_ - _VESSEL_FRONT_SIDE_THICKNESS_) / 2;
0164 #ifdef WITH_IRT2_SUPPORT
0165   double gvOffset = gas_volume_offset;
0166 #endif
0167   Tube gasTube(0.0, gas_volume_radius, gas_volume_length / 2);
0168   SubtractionSolid gasSolid(gasTube,
0169                             FlangeCut(description, gas_volume_length + 1 * mm, _FLANGE_CLEARANCE_));
0170   Volume gasVolume(detName + "_GasVol", gasSolid, vesselGas);
0171   pfRICH_volume.placeVolume(gasVolume, Position(0, 0, gas_volume_offset));
0172 #ifdef WITH_IRT2_SUPPORT
0173   {
0174     // FIXME: Z-location does not really matter here, right?;
0175     auto boundary =
0176         new FlatSurface(TVector3(0, 0, 0), sign * TVector3(1, 0, 0), TVector3(0, -1, 0));
0177
0178     auto radiator =
0179         geometry->SetContainerVolume(cdet, "GasVolume", 0, (G4LogicalVolume*)(0x10), 0, boundary);
0180     radiator->SetAlternativeMaterialName(vesselGasName.c_str());
0181   }
0182 #endif
0183
0184   auto _BUILDING_BLOCK_CLEARANCE_    = description.constant<double>("BUILDING_BLOCK_CLEARANCE");
0185   auto _VESSEL_INNER_WALL_THICKNESS_ = description.constant<double>("VESSEL_INNER_WALL_THICKNESS");
0186
0187   // To be used in boolean operations in several places;
0188   auto flange =
0189       FlangeCut(description, gas_volume_length + 1 * mm,
0190                 _FLANGE_CLEARANCE_ + _VESSEL_INNER_WALL_THICKNESS_ + _BUILDING_BLOCK_CLEARANCE_);
0191
0192   double gzOffset = -gas_volume_length / 2 + _BUILDING_BLOCK_CLEARANCE_;
0193
0194   auto _FLANGE_EPIPE_DIAMETER_ = description.constant<double>("FLANGE_EPIPE_DIAMETER");
0195   float m_r0min = _FLANGE_EPIPE_DIAMETER_ / 2 + _FLANGE_CLEARANCE_ + _VESSEL_INNER_WALL_THICKNESS_ +
0196                   _BUILDING_BLOCK_CLEARANCE_;
0197   float m_r0max = gas_volume_radius - _BUILDING_BLOCK_CLEARANCE_;
0198
0199   //
0200   // Aerogel
0201   //
0202   {
0203     auto _AEROGEL_INNER_WALL_THICKNESS_ =
0204         description.constant<double>("AEROGEL_INNER_WALL_THICKNESS");
0205     auto _AEROGEL_SEPARATOR_WALL_THICKNESS_ =
0206         description.constant<double>("AEROGEL_SEPARATOR_WALL_THICKNESS");
0207     auto _AEROGEL_OUTER_WALL_THICKNESS_ =
0208         description.constant<double>("AEROGEL_OUTER_WALL_THICKNESS");
0209
0210     auto aerogelMatName = aerogelElem.attr<std::string>(_Unicode(material));
0211     auto aerogelMat     = description.material(aerogelMatName);
0212     auto aerogelVis     = description.visAttributes(aerogelElem.attr<std::string>(_Unicode(vis)));
0213
0214     const int _AEROGEL_BAND_COUNT_            = 3;
0215     const unsigned adim[_AEROGEL_BAND_COUNT_] = {9, 14, 20};
0216     double rheight =
0217         (m_r0max - m_r0min - (_AEROGEL_BAND_COUNT_ - 1) * _AEROGEL_SEPARATOR_WALL_THICKNESS_ -
0218          _AEROGEL_INNER_WALL_THICKNESS_ - _AEROGEL_OUTER_WALL_THICKNESS_) /
0219         _AEROGEL_BAND_COUNT_;
0220
0221     for (unsigned ir = 0; ir < _AEROGEL_BAND_COUNT_; ir++) {
0222       double apitch     = 360 * degree / adim[ir];
0223       double aerogel_r0 = m_r0min + _AEROGEL_INNER_WALL_THICKNESS_ +
0224                           ir * (_AEROGEL_SEPARATOR_WALL_THICKNESS_ + rheight);
0225       double aerogel_r1 = aerogel_r0 + rheight;
0226       double rm         = (aerogel_r0 + aerogel_r1) / 2;
0227
0228       // Calculate angular space occupied by the spacers and by the tiles; no gas gaps for now;
0229       // assume that a wedge shape is good enough (GEANT visualization does not like boolean objects),
0230       // rather than creating constant thickness azimuthal spacers; just assume that spacer thickness is
0231       // _AEROGEL_FRAME_WALL_THICKNESS_ at r=rm;
0232       double l0   = 2 * M_PI * rm / adim[ir];
0233       double l1   = _AEROGEL_SEPARATOR_WALL_THICKNESS_;
0234       double lsum = l0 + l1;
0235
0236       // FIXME: names overlap in several places (?);
0237       double wd0 = (l0 / lsum) * (360 * degree / adim[ir]);
0238
0239       Tube agtube(aerogel_r0, aerogel_r1, aerogelThickness / 2, 0 * degree, wd0);
0240
0241       for (unsigned ia = 0; ia < adim[ir]; ia++) {
0242         Rotation3D r_aerogel_Z(RotationZYX(ia * apitch, 0.0, 0.0));
0243         Rotation3D r_aerogel_Zinv(RotationZYX(-1. * ia * apitch, 0.0, 0.0));
0244
0245         TString ag_name;
0246         ag_name.Form("PFRICH-aerogel-%d-%02d", ir, ia);
0247
0248         if (ir) {
0249           Volume agtubeVol(ag_name.Data(), agtube, aerogelMat);
0250           agtubeVol.setVisAttributes(aerogelVis);
0251           auto aerogelTilePlacement =
0252               Transform3D(r_aerogel_Z, Position(0.0, 0.0, gzOffset + aerogelThickness / 2));
0253           gasVolume.placeVolume(agtubeVol, aerogelTilePlacement);
0254         } else {
0255           Tube agtube_inner(aerogel_r0, aerogel_r1, aerogelThickness / 2, 0 * degree + ia * apitch,
0256                             wd0 + ia * apitch);
0257           SubtractionSolid agsub(agtube_inner, flange);
0258           Volume agsubtubeVol(ag_name.Data(), agsub, aerogelMat);
0259           agsubtubeVol.setVisAttributes(aerogelVis);
0260           auto aerogelTilePlacement = Transform3D(
0261               RotationZYX(0.0, 0.0, 0.0), Position(0.0, 0.0, gzOffset + aerogelThickness / 2));
0262           gasVolume.placeVolume(agsubtubeVol, aerogelTilePlacement);
0263         } //if
0264       } //for ia
0265     } // for ir
0266
0267 #ifdef WITH_IRT2_SUPPORT
0268     {
0269       TVector3 nx(1 * sign, 0, 0), ny(0, -1, 0);
0270
0271       auto surface = new FlatSurface(
0272           sign * (1 / mm) * TVector3(0, 0, fvOffset + gvOffset + gzOffset + aerogelThickness / 2),
0273           nx, ny);
0274
0275       auto radiator =
0276           geometry->AddFlatRadiator(cdet, "Aerogel", CherenkovDetector::Upstream, 0,
0277                                     (G4LogicalVolume*)(0x11), 0, surface, aerogelThickness / mm);
0278       radiator->SetAlternativeMaterialName(aerogelMatName.c_str());
0279     }
0280 #endif
0281
0282     // NB: there should be a small gap between aerogel and acrylic filter placed into the same
0283     // gas volume, otherwise IRT gets confused;
0284     gzOffset += aerogelThickness + _BUILDING_BLOCK_CLEARANCE_;
0285   }
0286
0287   //
0288   // Acrylic filter
0289   //
0290   {
0291     auto filterMatName = filterElem.attr<std::string>(_Unicode(material));
0292     auto filterMat     = description.material(filterMatName);
0293     auto filterVis     = description.visAttributes(filterElem.attr<std::string>(_Unicode(vis)));
0294
0295     Tube ac_tube(m_r0min, m_r0max, filterThickness / 2, 0 * degree, 360 * degree);
0296     SubtractionSolid ac_shape(ac_tube, flange);
0297     Volume acVol(detName + "-filter", ac_shape, filterMat);
0298     acVol.setVisAttributes(filterVis);
0299
0300     gasVolume.placeVolume(acVol, Position(0, 0, gzOffset + filterThickness / 2));
0301
0302 #ifdef WITH_IRT2_SUPPORT
0303     {
0304       TVector3 nx(1 * sign, 0, 0), ny(0, -1, 0);
0305
0306       auto surface = new FlatSurface(
0307           sign * (1 / mm) * TVector3(0, 0, fvOffset + gvOffset + gzOffset + filterThickness / 2),
0308           nx, ny);
0309
0310       auto radiator =
0311           geometry->AddFlatRadiator(cdet, "Acrylic", CherenkovDetector::Upstream, 0,
0312                                     (G4LogicalVolume*)(0x12), 0, surface, filterThickness / mm);
0313       radiator->SetAlternativeMaterialName(filterMatName.c_str());
0314     }
0315 #endif
0316   }
0317
0318 #ifdef WITH_IRT2_SUPPORT
0319   IRT2::OpticalBoundary* mboundaries[2] = {0, 0};
0320 #endif
0321
0322   //
0323   // Mirrors
0324   //
0325   {
0326     auto mirrorElem = detElem.child(_Unicode(mirror));
0327     auto mirrorMat  = description.material(mirrorElem.attr<std::string>(_Unicode(material)));
0328     auto mirrorVis  = description.visAttributes(mirrorElem.attr<std::string>(_Unicode(vis)));
0329     auto mirrorSurf = surfMgr.opticalSurface(mirrorElem.attr<std::string>(_Unicode(surface)));
0330
0331     auto _CONICAL_MIRROR_INNER_RADIUS_ =
0332         description.constant<double>("CONICAL_MIRROR_INNER_RADIUS");
0333     auto _CONICAL_MIRROR_OUTER_RADIUS_ =
0334         description.constant<double>("CONICAL_MIRROR_OUTER_RADIUS");
0335     auto _INNER_MIRROR_THICKNESS_ = description.constant<double>("INNER_MIRROR_THICKNESS");
0336     auto _OUTER_MIRROR_THICKNESS_ = description.constant<double>("OUTER_MIRROR_THICKNESS");
0337
0338     double mlen =
0339         gas_volume_length - 4 * _BUILDING_BLOCK_CLEARANCE_ - aerogelThickness - filterThickness;
0340     double mzoffset = (aerogelThickness + filterThickness + 2 * _BUILDING_BLOCK_CLEARANCE_) / 2;
0341
0342     double mirror_r0[2] = {m_r0min, m_r0max - _BUILDING_BLOCK_CLEARANCE_};
0343     double mirror_r1[2] = {_CONICAL_MIRROR_INNER_RADIUS_, _CONICAL_MIRROR_OUTER_RADIUS_};
0344
0345     for (unsigned im = 0; im < 2; im++) {
0346       double mirror_thickness = im ? _OUTER_MIRROR_THICKNESS_ : _INNER_MIRROR_THICKNESS_;
0347
0348       if (im) {
0349         Cone mirror_outer_cone_shape(mlen / 2.0, mirror_r0[im], mirror_r0[im] + mirror_thickness,
0350                                      mirror_r1[im], mirror_r1[im] + mirror_thickness);
0351
0352         Volume outer_mirrorVol(detName + "-outer-mirror", mirror_outer_cone_shape, mirrorMat);
0353         outer_mirrorVol.setVisAttributes(mirrorVis);
0354
0355         PlacedVolume mirror_outerPV =
0356             gasVolume.placeVolume(outer_mirrorVol, Position(0, 0, mzoffset));
0357         DetElement mirror_outerDE(sdet, "_outer_mirror_de", 0);
0358         mirror_outerDE.setPlacement(mirror_outerPV);
0359         SkinSurface mirrorSkin(description, mirror_outerDE, "outer_mirror_optical_surface_",
0360                                mirrorSurf, outer_mirrorVol);
0361         mirrorSkin.isValid();
0362
0363 #ifdef WITH_IRT2_SUPPORT
0364         auto msurface = new ConicalSurface(
0365             sign * (1 / mm) * TVector3(0, 0, fvOffset + gvOffset + mzoffset),
0366             sign * TVector3(0, 0, 1), mirror_r0[im] / mm, mirror_r1[im] / mm, mlen / mm);
0367
0368         mboundaries[im] =
0369             new IRT2::OpticalBoundary(cdet->GetRadiator("GasVolume"), msurface, false);
0370         // Need to store it in a separate call (?), see a comment in CherenkovDetector.h;
0371         cdet->StoreOpticalBoundary(mboundaries[im]);
0372 #endif
0373       } else {
0374
0375         Cone mirror_inner_cone_shape(mlen / 2., mirror_r0[im], mirror_r0[im] + mirror_thickness,
0376                                      mirror_r1[im], mirror_r1[im] + mirror_thickness);
0377
0378         SubtractionSolid mirror_inner_sub(mirror_inner_cone_shape, flange);
0379
0380         Volume inner_mirrorVol(detName + "-inner-mirror", mirror_inner_sub, mirrorMat);
0381         inner_mirrorVol.setVisAttributes(mirrorVis);
0382
0383         PlacedVolume mirror_innerPV =
0384             gasVolume.placeVolume(inner_mirrorVol, Position(0, 0, mzoffset));
0385         DetElement mirror_innerDE(sdet, "_inner_mirror_de", 0);
0386         mirror_innerDE.setPlacement(mirror_innerPV);
0387         SkinSurface mirrorSkin(description, mirror_innerDE, "inner_mirror_optical_surface_",
0388                                mirrorSurf, inner_mirrorVol);
0389         mirrorSkin.isValid();
0390
0391 #ifdef WITH_IRT2_SUPPORT
0392         auto msurface =
0393             new ConicalSurface(sign * (1 / mm) * TVector3(0, 0, fvOffset + gvOffset + mzoffset),
0394                                sign * TVector3(0, 0, 1), (mirror_r0[im] + mirror_thickness) / mm,
0395                                (mirror_r1[im] + mirror_thickness) / mm, mlen / mm);
0396         msurface->SetConvex();
0397
0398         mboundaries[im] =
0399             new IRT2::OpticalBoundary(cdet->GetRadiator("GasVolume"), msurface, false);
0400         // Need to store it in a separate call (?), see a comment in CherenkovDetector.h;
0401         cdet->StoreOpticalBoundary(mboundaries[im]);
0402 #endif
0403       }
0404     } //for im
0405   }
0406
0407   //
0408   // HRPPDs
0409   //
0410   {
0411     auto hrppdElem       = detElem.child(_Unicode(hrppd));
0412     auto HRPPD_WindowMat = description.material(hrppdElem.attr<std::string>(_Unicode(windowmat)));
0413     auto HRPPD_pcMat  = description.material(hrppdElem.attr<std::string>(_Unicode(photocathode)));
0414     auto HRPPD_MCPMat = description.material(hrppdElem.attr<std::string>(_Unicode(mcpmat)));
0415     auto HRPPD_PCBMat = description.material(hrppdElem.attr<std::string>(_Unicode(pcbmat)));
0416     auto HRPPD_PlatingMat = description.material(hrppdElem.attr<std::string>(_Unicode(plating)));
0417     auto HRPPD_CeramicMat = description.material(hrppdElem.attr<std::string>(_Unicode(ceramic)));
0418     auto pcVis   = description.visAttributes(hrppdElem.attr<std::string>(_Unicode(vis_pc)));
0419     auto wndVis  = description.visAttributes(hrppdElem.attr<std::string>(_Unicode(vis_window)));
0420     auto bodyVis = description.visAttributes(hrppdElem.attr<std::string>(_Unicode(vis_body)));
0421
0422     auto _HRPPD_CENTRAL_ROW_OFFSET_ = description.constant<double>("HRPPD_CENTRAL_ROW_OFFSET");
0423     auto _HRPPD_WINDOW_THICKNESS_   = description.constant<double>("HRPPD_WINDOW_THICKNESS");
0424     auto _HRPPD_CONTAINER_VOLUME_HEIGHT_ =
0425         description.constant<double>("HRPPD_CONTAINER_VOLUME_HEIGHT");
0426     auto _HRPPD_INSTALLATION_GAP_ = description.constant<double>("HRPPD_INSTALLATION_GAP");
0427
0428     auto _HRPPD_TILE_SIZE_        = description.constant<double>("HRPPD_TILE_SIZE");
0429     auto _HRPPD_OPEN_AREA_SIZE_   = description.constant<double>("HRPPD_OPEN_AREA_SIZE");
0430     auto _HRPPD_ACTIVE_AREA_SIZE_ = description.constant<double>("HRPPD_ACTIVE_AREA_SIZE");
0431     auto _HRPPD_CERAMIC_BODY_THICKNESS_ =
0432         description.constant<double>("HRPPD_CERAMIC_BODY_THICKNESS");
0433     auto _HRPPD_PHOTOCATHODE_THICKNESS_ =
0434         description.constant<double>("HRPPD_PHOTOCATHODE_THICKNESS");
0435     auto _HRPPD_BASEPLATE_THICKNESS_ = description.constant<double>("HRPPD_BASEPLATE_THICKNESS");
0436     auto _HRPPD_PLATING_LAYER_THICKNESS_ =
0437         description.constant<double>("HRPPD_PLATING_LAYER_THICKNESS");
0438     auto _EFFECTIVE_MCP_THICKNESS_ = description.constant<double>("EFFECTIVE_MCP_THICKNESS");
0439 #ifdef WITH_IRT2_SUPPORT
0440     auto _HRPPD_COLLECTION_EFFICIENCY_ =
0441         description.constant<double>("HRPPD_COLLECTION_EFFICIENCY");
0442 #endif
0443
0444 #ifdef WITH_IRT2_SUPPORT
0445     // [0,0]: have neither access to G4VSolid nor to G4Material; IRT code does not care; fine;
0446     auto pd = new IRT2::CherenkovPhotonDetector(0, 0);
0447
0448     // FIXME: '0' stands for the unknown (and irrelevant) G4LogicalVolume;
0449     geometry->AddPhotonDetector(cdet, (G4LogicalVolume*)0x10, pd);
0450
0451     // Cannot access GEANT shapes in the reconstruction code -> store this value;
0452     pd->SetActiveAreaSize(_HRPPD_ACTIVE_AREA_SIZE_ / mm);
0453
0454     pd->SetGeometricEfficiency(_HRPPD_COLLECTION_EFFICIENCY_);
0455 #endif
0456
0457 #ifdef WITH_IRT2_SUPPORT
0458     {
0459       TVector3 nx(1 * sign, 0, 0), ny(0, -1, 0);
0460
0461       // For now assume it is a unique surface, same for all HRPPDs;
0462       auto surface =
0463           new FlatSurface(sign * (1 / mm) *
0464                               TVector3(0, 0,
0465                                        fvOffset + _FIDUCIAL_VOLUME_LENGTH_ / 2 -
0466                                            _SENSOR_AREA_LENGTH_ + _HRPPD_WINDOW_THICKNESS_ / 2),
0467                           nx, ny);
0468
0469       auto radiator = geometry->AddFlatRadiator(cdet, "QuartzWindow", CherenkovDetector::Downstream,
0470                                                 0, (G4LogicalVolume*)(0x13), 0, surface,
0471                                                 _HRPPD_WINDOW_THICKNESS_ / mm);
0472       radiator->SetAlternativeMaterialName("AirOptical");
0473     }
0474 #endif
0475
0476     // Create a vector of HRPPD XY-coordinates in a separate loop;
0477     std::vector<std::pair<TVector2, bool>> coord;
0478     {
0479       unsigned const hdim              = 9;
0480       const unsigned flags[hdim][hdim] = {
0481           // NB: WYSIWIG fashion; well, it is top/bottom and left/right symmetric;
0482           // clang-format off
0483         {0, 0, 1, 1, 1, 1, 1, 0, 0},
0484         {0, 1, 1, 1, 1, 1, 1, 1, 0},
0485         {1, 1, 1, 1, 1, 1, 1, 1, 1},
0486         {1, 1, 1, 1, 2, 1, 1, 1, 1},
0487         {3, 3, 3, 4, 0, 2, 1, 1, 1},
0488         {1, 1, 1, 1, 2, 1, 1, 1, 1},
0489         {1, 1, 1, 1, 1, 1, 1, 1, 1},
0490         {0, 1, 1, 1, 1, 1, 1, 1, 0},
0491         {0, 0, 1, 1, 1, 1, 1, 0, 0}};
0492       // clang-format on
0493
0494       for (unsigned ix = 0; ix < hdim; ix++) {
0495         double xOffset = (_HRPPD_TILE_SIZE_ + _HRPPD_INSTALLATION_GAP_) * (ix - (hdim - 1) / 2.);
0496
0497         for (unsigned iy = 0; iy < hdim; iy++) {
0498           double yOffset = (_HRPPD_TILE_SIZE_ + _HRPPD_INSTALLATION_GAP_) * (iy - (hdim - 1) / 2.);
0499           unsigned flag  = flags[hdim - iy - 1][ix];
0500
0501           if (!flag)
0502             continue;
0503
0504           double qxOffset = xOffset + (flag >= 3 ? -_HRPPD_CENTRAL_ROW_OFFSET_ : 0.0);
0505           coord.push_back(std::make_pair(TVector2(qxOffset, yOffset), flag % 2));
0506         } //for iy
0507       } //for ix
0508     }
0509
0510     unsigned imod = 0;
0511
0512     //
0513     // It looks like each HRPPD should be a new object rather than a copy of the same volume,
0514     // because otherwise one cannot make photocathodes sensitive (this is done on a PlacedVolume
0515     // level);
0516     //
0517     for (auto xyptr : coord) {
0518       auto& xy = xyptr.first;
0519
0520       uint64_t sensorID = 0x0;
0521
0522       //
0523       // HRPPD container volume
0524       //
0525       Box hrppd_Solid(_HRPPD_TILE_SIZE_ / 2, _HRPPD_TILE_SIZE_ / 2,
0526                       _HRPPD_CONTAINER_VOLUME_HEIGHT_ / 2);
0527       TString hrppdName;
0528       hrppdName.Form("%s-hrppd-%02d", detName.c_str(), imod);
0529       // FIXME: may want to use AirOptical here, but then return a dummy absorber
0530       // layer behind the actual photocathode;
0531       Volume hrppdVol_air(hrppdName.Data(), hrppd_Solid, air);
0532
0533       // A running variable to pack layers one after the other one;
0534       double accu = -_HRPPD_CONTAINER_VOLUME_HEIGHT_ / 2;
0535
0536       //
0537       // Quartz Window
0538       //
0539       Box wnd_Solid(_HRPPD_TILE_SIZE_ / 2, _HRPPD_TILE_SIZE_ / 2, _HRPPD_WINDOW_THICKNESS_ / 2);
0540       TString wndName;
0541       wndName.Form("%s-window-%02d", detName.c_str(), imod);
0542       Volume wndVol(wndName.Data(), wnd_Solid, HRPPD_WindowMat);
0543       wndVol.setVisAttributes(wndVis);
0544       hrppdVol_air.placeVolume(wndVol, Position(0, 0, accu + _HRPPD_WINDOW_THICKNESS_ / 2));
0545
0546       accu += _HRPPD_WINDOW_THICKNESS_;
0547
0548       //
0549       // Photocathode layer (sensitive volume)
0550       //
0551       {
0552         auto pcBox = Box(_HRPPD_ACTIVE_AREA_SIZE_ / 2, _HRPPD_ACTIVE_AREA_SIZE_ / 2,
0553                          _HRPPD_PHOTOCATHODE_THICKNESS_ / 2);
0554         TString pcName;
0555         pcName.Form("%s-photocathode-%02d", detName.c_str(), imod);
0556         Volume pcVol(pcName.Data(), pcBox, HRPPD_pcMat);
0557         pcVol.setSensitiveDetector(sens);
0558
0559         pcVol.setVisAttributes(pcVis);
0560         PlacedVolume pcPV = hrppdVol_air.placeVolume(
0561             pcVol, Position(0.0, 0.0, accu + _HRPPD_PHOTOCATHODE_THICKNESS_ / 2));
0562         {
0563           pcPV.addPhysVolID("hrppd", imod);
0564
0565           // sensor DetElement
0566           sensorID = encodeSensorID(pcPV.volIDs());
0567           TString deName;
0568           deName.Form("%s-sensor-%02d", detName.c_str(), imod);
0569           DetElement pcDE(sdet, deName.Data(), sensorID);
0570           pcDE.setPlacement(pcPV);
0571         }
0572
0573         //
0574         // A fake absorber layer behind the photocathode; FIXME: make sure that reflection
0575         // on the window and photocathode boundary still works correctly (no fake volume as
0576         // in a standalone code);
0577         //
0578         {
0579           TString absName;
0580           absName.Form("%s-absorber-%02d", detName.c_str(), imod);
0581           // Recycle the same pcBox shape; do not mind to use PCB material;
0582           Volume absVol(absName.Data(), pcBox, HRPPD_PCBMat);
0583           hrppdVol_air.placeVolume(absVol, Position(0.0, 0.0,
0584                                                     accu + _HRPPD_PHOTOCATHODE_THICKNESS_ +
0585                                                         _HRPPD_PHOTOCATHODE_THICKNESS_ / 2));
0586         }
0587       }
0588
0589       //
0590       // Ceramic body (sidewall and anode)
0591       //
0592       {
0593         Box cerbox(_HRPPD_TILE_SIZE_ / 2, _HRPPD_TILE_SIZE_ / 2,
0594                    _HRPPD_CERAMIC_BODY_THICKNESS_ / 2);
0595         Box cutbox(_HRPPD_OPEN_AREA_SIZE_ / 2, _HRPPD_OPEN_AREA_SIZE_ / 2,
0596                    _HRPPD_CERAMIC_BODY_THICKNESS_ / 2);
0597
0598         SubtractionSolid ceramic(cerbox, cutbox, Position(0, 0, -_HRPPD_BASEPLATE_THICKNESS_));
0599
0600         TString cerName;
0601         cerName.Form("%s-ceramic-%02d", detName.c_str(), imod);
0602         Volume ceramicVol(cerName.Data(), ceramic, HRPPD_CeramicMat);
0603
0604         ceramicVol.setVisAttributes(bodyVis);
0605         hrppdVol_air.placeVolume(ceramicVol,
0606                                  Position(0.0, 0.0, accu + _HRPPD_CERAMIC_BODY_THICKNESS_ / 2));
0607       }
0608
0609       //
0610       // Effective anode plating layer
0611       //
0612       {
0613         Box plating_solid(_HRPPD_OPEN_AREA_SIZE_ / 2, _HRPPD_OPEN_AREA_SIZE_ / 2,
0614                           _HRPPD_PLATING_LAYER_THICKNESS_ / 2);
0615         TString pltName;
0616         pltName.Form("%s-plating-%02d", detName.c_str(), imod);
0617         Volume platingVol(pltName.Data(), plating_solid, HRPPD_PlatingMat);
0618         // Place somewhere in the middle of the ceramic body gap;
0619         hrppdVol_air.placeVolume(platingVol,
0620                                  Position(0.0, 0.0, accu + _HRPPD_CERAMIC_BODY_THICKNESS_ / 2));
0621       }
0622
0623       //
0624       // Effective MCP layer
0625       //
0626       {
0627         Box mcp_solid(_HRPPD_OPEN_AREA_SIZE_ / 2, _HRPPD_OPEN_AREA_SIZE_ / 2,
0628                       _EFFECTIVE_MCP_THICKNESS_ / 2);
0629         TString mcpName;
0630         mcpName.Form("%s-mcp-%02d", detName.c_str(), imod);
0631         Volume mcpVol(mcpName.Data(), mcp_solid, HRPPD_MCPMat);
0632         hrppdVol_air.placeVolume(mcpVol, Position(0.0, 0.0,
0633                                                   accu + _HRPPD_CERAMIC_BODY_THICKNESS_ / 2 +
0634                                                       _HRPPD_PLATING_LAYER_THICKNESS_ +
0635                                                       _EFFECTIVE_MCP_THICKNESS_ / 2));
0636       }
0637
0638       accu += _HRPPD_CERAMIC_BODY_THICKNESS_;
0639
0640       //
0641       // PCB
0642       //
0643       {
0644         auto _READOUT_PCB_THICKNESS_ = description.constant<double>("READOUT_PCB_THICKNESS");
0645         auto _READOUT_PCB_SIZE_      = description.constant<double>("READOUT_PCB_SIZE");
0646
0647         Box pcb_solid(_READOUT_PCB_SIZE_ / 2, _READOUT_PCB_SIZE_ / 2, _READOUT_PCB_THICKNESS_ / 2);
0648         TString pcbName;
0649         pcbName.Form("%s-pcb-%02d", detName.c_str(), imod);
0650         Volume pcbVol(pcbName.Data(), pcb_solid, HRPPD_PCBMat);
0651         hrppdVol_air.placeVolume(pcbVol, Position(0.0, 0.0, accu + _READOUT_PCB_THICKNESS_ / 2));
0652       }
0653
0654       // Eventually place the whole HRPPD container volume;
0655       double dz =
0656           _FIDUCIAL_VOLUME_LENGTH_ / 2 - _SENSOR_AREA_LENGTH_ + _HRPPD_CONTAINER_VOLUME_HEIGHT_ / 2;
0657       pfRICH_volume.placeVolume(hrppdVol_air, Position(xy.X(), xy.Y(), dz));
0658
0659 #ifdef WITH_IRT2_SUPPORT
0660       {
0661         // Photocathode surface;
0662         double xOffset = xy.X(), yOffset = xy.Y();
0663         auto surface = new FlatSurface(
0664             (1 / mm) *
0665                 TVector3(sign * xOffset, yOffset,
0666                          sign * (fvOffset + _FIDUCIAL_VOLUME_LENGTH_ / 2 - _SENSOR_AREA_LENGTH_ +
0667                                  _HRPPD_WINDOW_THICKNESS_ + _HRPPD_PHOTOCATHODE_THICKNESS_ / 2)),
0668             TVector3(1 * sign, 0, 0), TVector3(0, -1, 0));
0669
0670         {
0671           // '0': pfRICH has no division in sectors (unlike e.g. dRICH);
0672           unsigned sector = 0;
0673
0674           for (unsigned iq = 0; iq < 4; iq++) {
0675             auto irt = pd->AllocateIRT(sector, sensorID);
0676
0677             // Aerogel and acrylic;
0678             if (cdet->m_OpticalBoundaries[CherenkovDetector::Upstream].find(sector) !=
0679                 cdet->m_OpticalBoundaries[CherenkovDetector::Upstream].end())
0680               for (auto boundary : cdet->m_OpticalBoundaries[CherenkovDetector::Upstream][sector])
0681                 irt->AddOpticalBoundary(boundary);
0682
0683             switch (iq) {
0684             case 0:
0685               // Direct hit;
0686               break;
0687             case 1:
0688             case 2:
0689               // Reflection on either inner or outer mirrors;
0690               irt->AddOpticalBoundary(mboundaries[iq - 1]);
0691               break;
0692             case 3:
0693               // Reflection on outer, then on inner mirror; happens at large angles; if the pyramids are
0694               // too high, these photons will undergo more reflections, and cannot be saved;
0695               irt->AddOpticalBoundary(mboundaries[1]);
0696               irt->AddOpticalBoundary(mboundaries[0]);
0697               break;
0698             } //switch
0699
0700             // Fused silica windows;
0701             if (cdet->m_OpticalBoundaries[CherenkovDetector::Downstream].find(sector) !=
0702                 cdet->m_OpticalBoundaries[CherenkovDetector::Downstream].end())
0703               for (auto boundary : cdet->m_OpticalBoundaries[CherenkovDetector::Downstream][sector])
0704                 irt->AddOpticalBoundary(boundary);
0705
0706             // Terminate the optical path;
0707             pd->AddItselfToOpticalBoundaries(irt, surface);
0708           } //for iq
0709         }
0710       }
0711 #endif
0712
0713       imod++;
0714     } //for coord
0715   }
0716
0717   return sdet;
0718 } // createDetector()
0719
0720 // -------------------------------------------------------------------------------------
0721
0722 // clang-format off
0723 DECLARE_DETELEMENT(epic_PFRICH, createDetector)