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
0139   Volume mother                 = description.pickMotherVolume(sdet);
0140   auto _FIDUCIAL_VOLUME_OFFSET_ = description.constant<double>("FIDUCIAL_VOLUME_OFFSET");
0141   Transform3D transform(RotationZYX(0, M_PI, 0), Position(0, 0, _FIDUCIAL_VOLUME_OFFSET_));
0142   PlacedVolume pv = mother.placeVolume(pfRICH_volume, transform);
0143   if (id != 0)
0144     pv.addPhysVolID("system", id);
0145   sdet.setPlacement(pv);
0146
0147   //
0148   // Gas volume
0149   //
0150   auto _VESSEL_FRONT_SIDE_THICKNESS_ = description.constant<double>("VESSEL_FRONT_SIDE_THICKNESS");
0151   auto _SENSOR_AREA_LENGTH_          = description.constant<double>("SENSOR_AREA_LENGTH");
0152   auto _VESSEL_OUTER_WALL_THICKNESS_ = description.constant<double>("VESSEL_OUTER_WALL_THICKNESS");
0153
0154   // FIXME: do it better later;
0155 #ifdef WITH_IRT2_SUPPORT
0156   double fvOffset = fabs(_FIDUCIAL_VOLUME_OFFSET_);
0157 #endif
0158
0159   double gas_volume_length =
0160       _FIDUCIAL_VOLUME_LENGTH_ - _VESSEL_FRONT_SIDE_THICKNESS_ - _SENSOR_AREA_LENGTH_;
0161   double gas_volume_radius = _VESSEL_OUTER_RADIUS_ - _VESSEL_OUTER_WALL_THICKNESS_;
0162   double gas_volume_offset = -(_SENSOR_AREA_LENGTH_ - _VESSEL_FRONT_SIDE_THICKNESS_) / 2;
0163 #ifdef WITH_IRT2_SUPPORT
0164   double gvOffset = gas_volume_offset;
0165 #endif
0166   Tube gasTube(0.0, gas_volume_radius, gas_volume_length / 2);
0167   SubtractionSolid gasSolid(gasTube,
0168                             FlangeCut(description, gas_volume_length + 1 * mm, _FLANGE_CLEARANCE_));
0169   Volume gasVolume(detName + "_GasVol", gasSolid, vesselGas);
0170   pfRICH_volume.placeVolume(gasVolume, Position(0, 0, gas_volume_offset));
0171 #ifdef WITH_IRT2_SUPPORT
0172   {
0173     // FIXME: Z-location does not really matter here, right?;
0174     auto boundary =
0175         new FlatSurface(TVector3(0, 0, 0), sign * TVector3(1, 0, 0), TVector3(0, -1, 0));
0176
0177     auto radiator =
0178         geometry->SetContainerVolume(cdet, "GasVolume", 0, (G4LogicalVolume*)(0x0), 0, boundary);
0179     radiator->SetAlternativeMaterialName(vesselGasName.c_str());
0180   }
0181 #endif
0182
0183   auto _BUILDING_BLOCK_CLEARANCE_    = description.constant<double>("BUILDING_BLOCK_CLEARANCE");
0184   auto _VESSEL_INNER_WALL_THICKNESS_ = description.constant<double>("VESSEL_INNER_WALL_THICKNESS");
0185
0186   // To be used in boolean operations in several places;
0187   auto flange =
0188       FlangeCut(description, gas_volume_length + 1 * mm,
0189                 _FLANGE_CLEARANCE_ + _VESSEL_INNER_WALL_THICKNESS_ + _BUILDING_BLOCK_CLEARANCE_);
0190
0191   double gzOffset = -gas_volume_length / 2 + _BUILDING_BLOCK_CLEARANCE_;
0192
0193   auto _FLANGE_EPIPE_DIAMETER_ = description.constant<double>("FLANGE_EPIPE_DIAMETER");
0194   float m_r0min = _FLANGE_EPIPE_DIAMETER_ / 2 + _FLANGE_CLEARANCE_ + _VESSEL_INNER_WALL_THICKNESS_ +
0195                   _BUILDING_BLOCK_CLEARANCE_;
0196   float m_r0max = gas_volume_radius - _BUILDING_BLOCK_CLEARANCE_;
0197
0198   //
0199   // Aerogel
0200   //
0201   {
0202     auto _AEROGEL_INNER_WALL_THICKNESS_ =
0203         description.constant<double>("AEROGEL_INNER_WALL_THICKNESS");
0204     auto _AEROGEL_SEPARATOR_WALL_THICKNESS_ =
0205         description.constant<double>("AEROGEL_SEPARATOR_WALL_THICKNESS");
0206     auto _AEROGEL_OUTER_WALL_THICKNESS_ =
0207         description.constant<double>("AEROGEL_OUTER_WALL_THICKNESS");
0208
0209     auto aerogelMatName = aerogelElem.attr<std::string>(_Unicode(material));
0210     auto aerogelMat     = description.material(aerogelMatName);
0211     auto aerogelVis     = description.visAttributes(aerogelElem.attr<std::string>(_Unicode(vis)));
0212
0213     const int _AEROGEL_BAND_COUNT_            = 3;
0214     const unsigned adim[_AEROGEL_BAND_COUNT_] = {9, 14, 20};
0215     double rheight =
0216         (m_r0max - m_r0min - (_AEROGEL_BAND_COUNT_ - 1) * _AEROGEL_SEPARATOR_WALL_THICKNESS_ -
0217          _AEROGEL_INNER_WALL_THICKNESS_ - _AEROGEL_OUTER_WALL_THICKNESS_) /
0218         _AEROGEL_BAND_COUNT_;
0219
0220     for (unsigned ir = 0; ir < _AEROGEL_BAND_COUNT_; ir++) {
0221       double apitch     = 360 * degree / adim[ir];
0222       double aerogel_r0 = m_r0min + _AEROGEL_INNER_WALL_THICKNESS_ +
0223                           ir * (_AEROGEL_SEPARATOR_WALL_THICKNESS_ + rheight);
0224       double aerogel_r1 = aerogel_r0 + rheight;
0225       double rm         = (aerogel_r0 + aerogel_r1) / 2;
0226
0227       // Calculate angular space occupied by the spacers and by the tiles; no gas gaps for now;
0228       // assume that a wedge shape is good enough (GEANT visualization does not like boolean objects),
0229       // rather than creating constant thickness azimuthal spacers; just assume that spacer thickness is
0230       // _AEROGEL_FRAME_WALL_THICKNESS_ at r=rm;
0231       double l0   = 2 * M_PI * rm / adim[ir];
0232       double l1   = _AEROGEL_SEPARATOR_WALL_THICKNESS_;
0233       double lsum = l0 + l1;
0234
0235       // FIXME: names overlap in several places (?);
0236       double wd0 = (l0 / lsum) * (360 * degree / adim[ir]);
0237
0238       Tube agtube(aerogel_r0, aerogel_r1, aerogelThickness / 2, 0 * degree, wd0);
0239
0240       for (unsigned ia = 0; ia < adim[ir]; ia++) {
0241         Rotation3D r_aerogel_Z(RotationZYX(ia * apitch, 0.0, 0.0));
0242         Rotation3D r_aerogel_Zinv(RotationZYX(-1. * ia * apitch, 0.0, 0.0));
0243
0244         TString ag_name;
0245         ag_name.Form("PFRICH-aerogel-%d-%02d", ir, ia);
0246
0247         if (ir) {
0248           Volume agtubeVol(ag_name.Data(), agtube, aerogelMat);
0249           agtubeVol.setVisAttributes(aerogelVis);
0250           auto aerogelTilePlacement =
0251               Transform3D(r_aerogel_Z, Position(0.0, 0.0, gzOffset + aerogelThickness / 2));
0252           gasVolume.placeVolume(agtubeVol, aerogelTilePlacement);
0253         } else {
0254           Tube agtube_inner(aerogel_r0, aerogel_r1, aerogelThickness / 2, 0 * degree + ia * apitch,
0255                             wd0 + ia * apitch);
0256           SubtractionSolid agsub(agtube_inner, flange);
0257           Volume agsubtubeVol(ag_name.Data(), agsub, aerogelMat);
0258           agsubtubeVol.setVisAttributes(aerogelVis);
0259           auto aerogelTilePlacement = Transform3D(
0260               RotationZYX(0.0, 0.0, 0.0), Position(0.0, 0.0, gzOffset + aerogelThickness / 2));
0261           gasVolume.placeVolume(agsubtubeVol, aerogelTilePlacement);
0262         } //if
0263       } //for ia
0264     } // for ir
0265
0266 #ifdef WITH_IRT2_SUPPORT
0267     {
0268       TVector3 nx(1 * sign, 0, 0), ny(0, -1, 0);
0269
0270       auto surface = new FlatSurface(
0271           sign * (1 / mm) * TVector3(0, 0, fvOffset + gvOffset + gzOffset + aerogelThickness / 2),
0272           nx, ny);
0273
0274       auto radiator =
0275           geometry->AddFlatRadiator(cdet, "Aerogel", CherenkovDetector::Upstream, 0,
0276                                     (G4LogicalVolume*)(0x1), 0, surface, aerogelThickness / mm);
0277       radiator->SetAlternativeMaterialName(aerogelMatName.c_str());
0278     }
0279 #endif
0280
0281     // NB: there should be a small gap between aerogel and acrylic filter placed into the same
0282     // gas volume, otherwise IRT gets confused;
0283     gzOffset += aerogelThickness + _BUILDING_BLOCK_CLEARANCE_;
0284   }
0285
0286   //
0287   // Acrylic filter
0288   //
0289   {
0290     auto filterMatName = filterElem.attr<std::string>(_Unicode(material));
0291     auto filterMat     = description.material(filterMatName);
0292     auto filterVis     = description.visAttributes(filterElem.attr<std::string>(_Unicode(vis)));
0293
0294     Tube ac_tube(m_r0min, m_r0max, filterThickness / 2, 0 * degree, 360 * degree);
0295     SubtractionSolid ac_shape(ac_tube, flange);
0296     Volume acVol(detName + "-filter", ac_shape, filterMat);
0297     acVol.setVisAttributes(filterVis);
0298
0299     gasVolume.placeVolume(acVol, Position(0, 0, gzOffset + filterThickness / 2));
0300
0301 #ifdef WITH_IRT2_SUPPORT
0302     {
0303       TVector3 nx(1 * sign, 0, 0), ny(0, -1, 0);
0304
0305       auto surface = new FlatSurface(
0306           sign * (1 / mm) * TVector3(0, 0, fvOffset + gvOffset + gzOffset + filterThickness / 2),
0307           nx, ny);
0308
0309       auto radiator =
0310           geometry->AddFlatRadiator(cdet, "Acrylic", CherenkovDetector::Upstream, 0,
0311                                     (G4LogicalVolume*)(0x2), 0, surface, filterThickness / mm);
0312       radiator->SetAlternativeMaterialName(filterMatName.c_str());
0313     }
0314 #endif
0315   }
0316
0317 #ifdef WITH_IRT2_SUPPORT
0318   IRT2::OpticalBoundary* mboundaries[2] = {0, 0};
0319 #endif
0320
0321   //
0322   // Mirrors
0323   //
0324   {
0325     auto mirrorElem = detElem.child(_Unicode(mirror));
0326     auto mirrorMat  = description.material(mirrorElem.attr<std::string>(_Unicode(material)));
0327     auto mirrorVis  = description.visAttributes(mirrorElem.attr<std::string>(_Unicode(vis)));
0328     auto mirrorSurf = surfMgr.opticalSurface(mirrorElem.attr<std::string>(_Unicode(surface)));
0329
0330     auto _CONICAL_MIRROR_INNER_RADIUS_ =
0331         description.constant<double>("CONICAL_MIRROR_INNER_RADIUS");
0332     auto _CONICAL_MIRROR_OUTER_RADIUS_ =
0333         description.constant<double>("CONICAL_MIRROR_OUTER_RADIUS");
0334     auto _INNER_MIRROR_THICKNESS_ = description.constant<double>("INNER_MIRROR_THICKNESS");
0335     auto _OUTER_MIRROR_THICKNESS_ = description.constant<double>("OUTER_MIRROR_THICKNESS");
0336
0337     double mlen =
0338         gas_volume_length - 4 * _BUILDING_BLOCK_CLEARANCE_ - aerogelThickness - filterThickness;
0339     double mzoffset = (aerogelThickness + filterThickness + 2 * _BUILDING_BLOCK_CLEARANCE_) / 2;
0340
0341     double mirror_r0[2] = {m_r0min, m_r0max - _BUILDING_BLOCK_CLEARANCE_};
0342     double mirror_r1[2] = {_CONICAL_MIRROR_INNER_RADIUS_, _CONICAL_MIRROR_OUTER_RADIUS_};
0343
0344     for (unsigned im = 0; im < 2; im++) {
0345       double mirror_thickness = im ? _OUTER_MIRROR_THICKNESS_ : _INNER_MIRROR_THICKNESS_;
0346
0347       if (im) {
0348         Cone mirror_outer_cone_shape(mlen / 2.0, mirror_r0[im], mirror_r0[im] + mirror_thickness,
0349                                      mirror_r1[im], mirror_r1[im] + mirror_thickness);
0350
0351         Volume outer_mirrorVol(detName + "-outer-mirror", mirror_outer_cone_shape, mirrorMat);
0352         outer_mirrorVol.setVisAttributes(mirrorVis);
0353
0354         PlacedVolume mirror_outerPV =
0355             gasVolume.placeVolume(outer_mirrorVol, Position(0, 0, mzoffset));
0356         DetElement mirror_outerDE(sdet, "_outer_mirror_de", 0);
0357         mirror_outerDE.setPlacement(mirror_outerPV);
0358         SkinSurface mirrorSkin(description, mirror_outerDE, "outer_mirror_optical_surface_",
0359                                mirrorSurf, outer_mirrorVol);
0360         mirrorSkin.isValid();
0361
0362 #ifdef WITH_IRT2_SUPPORT
0363         auto msurface = new ConicalSurface(
0364             sign * (1 / mm) * TVector3(0, 0, fvOffset + gvOffset + mzoffset),
0365             sign * TVector3(0, 0, 1), mirror_r0[im] / mm, mirror_r1[im] / mm, mlen / mm);
0366
0367         mboundaries[im] =
0368             new IRT2::OpticalBoundary(cdet->GetRadiator("GasVolume"), msurface, false);
0369         // Need to store it in a separate call (?), see a comment in CherenkovDetector.h;
0370         cdet->StoreOpticalBoundary(mboundaries[im]);
0371 #endif
0372       } else {
0373
0374         Cone mirror_inner_cone_shape(mlen / 2., mirror_r0[im], mirror_r0[im] + mirror_thickness,
0375                                      mirror_r1[im], mirror_r1[im] + mirror_thickness);
0376
0377         SubtractionSolid mirror_inner_sub(mirror_inner_cone_shape, flange);
0378
0379         Volume inner_mirrorVol(detName + "-inner-mirror", mirror_inner_sub, mirrorMat);
0380         inner_mirrorVol.setVisAttributes(mirrorVis);
0381
0382         PlacedVolume mirror_innerPV =
0383             gasVolume.placeVolume(inner_mirrorVol, Position(0, 0, mzoffset));
0384         DetElement mirror_innerDE(sdet, "_inner_mirror_de", 0);
0385         mirror_innerDE.setPlacement(mirror_innerPV);
0386         SkinSurface mirrorSkin(description, mirror_innerDE, "inner_mirror_optical_surface_",
0387                                mirrorSurf, inner_mirrorVol);
0388         mirrorSkin.isValid();
0389
0390 #ifdef WITH_IRT2_SUPPORT
0391         auto msurface =
0392             new ConicalSurface(sign * (1 / mm) * TVector3(0, 0, fvOffset + gvOffset + mzoffset),
0393                                sign * TVector3(0, 0, 1), (mirror_r0[im] + mirror_thickness) / mm,
0394                                (mirror_r1[im] + mirror_thickness) / mm, mlen / mm);
0395         msurface->SetConvex();
0396
0397         mboundaries[im] =
0398             new IRT2::OpticalBoundary(cdet->GetRadiator("GasVolume"), msurface, false);
0399         // Need to store it in a separate call (?), see a comment in CherenkovDetector.h;
0400         cdet->StoreOpticalBoundary(mboundaries[im]);
0401 #endif
0402       }
0403     } //for im
0404   }
0405
0406   //
0407   // HRPPDs
0408   //
0409   {
0410     auto hrppdElem       = detElem.child(_Unicode(hrppd));
0411     auto HRPPD_WindowMat = description.material(hrppdElem.attr<std::string>(_Unicode(windowmat)));
0412     auto HRPPD_pcMat  = description.material(hrppdElem.attr<std::string>(_Unicode(photocathode)));
0413     auto HRPPD_MCPMat = description.material(hrppdElem.attr<std::string>(_Unicode(mcpmat)));
0414     auto HRPPD_PCBMat = description.material(hrppdElem.attr<std::string>(_Unicode(pcbmat)));
0415     auto HRPPD_PlatingMat = description.material(hrppdElem.attr<std::string>(_Unicode(plating)));
0416     auto HRPPD_CeramicMat = description.material(hrppdElem.attr<std::string>(_Unicode(ceramic)));
0417     auto pcVis   = description.visAttributes(hrppdElem.attr<std::string>(_Unicode(vis_pc)));
0418     auto wndVis  = description.visAttributes(hrppdElem.attr<std::string>(_Unicode(vis_window)));
0419     auto bodyVis = description.visAttributes(hrppdElem.attr<std::string>(_Unicode(vis_body)));
0420
0421     auto _HRPPD_CENTRAL_ROW_OFFSET_ = description.constant<double>("HRPPD_CENTRAL_ROW_OFFSET");
0422     auto _HRPPD_WINDOW_THICKNESS_   = description.constant<double>("HRPPD_WINDOW_THICKNESS");
0423     auto _HRPPD_CONTAINER_VOLUME_HEIGHT_ =
0424         description.constant<double>("HRPPD_CONTAINER_VOLUME_HEIGHT");
0425     auto _HRPPD_INSTALLATION_GAP_ = description.constant<double>("HRPPD_INSTALLATION_GAP");
0426
0427     auto _HRPPD_TILE_SIZE_        = description.constant<double>("HRPPD_TILE_SIZE");
0428     auto _HRPPD_OPEN_AREA_SIZE_   = description.constant<double>("HRPPD_OPEN_AREA_SIZE");
0429     auto _HRPPD_ACTIVE_AREA_SIZE_ = description.constant<double>("HRPPD_ACTIVE_AREA_SIZE");
0430     auto _HRPPD_CERAMIC_BODY_THICKNESS_ =
0431         description.constant<double>("HRPPD_CERAMIC_BODY_THICKNESS");
0432     auto _HRPPD_PHOTOCATHODE_THICKNESS_ =
0433         description.constant<double>("HRPPD_PHOTOCATHODE_THICKNESS");
0434     auto _HRPPD_BASEPLATE_THICKNESS_ = description.constant<double>("HRPPD_BASEPLATE_THICKNESS");
0435     auto _HRPPD_PLATING_LAYER_THICKNESS_ =
0436         description.constant<double>("HRPPD_PLATING_LAYER_THICKNESS");
0437     auto _EFFECTIVE_MCP_THICKNESS_ = description.constant<double>("EFFECTIVE_MCP_THICKNESS");
0438 #ifdef WITH_IRT2_SUPPORT
0439     auto _HRPPD_COLLECTION_EFFICIENCY_ =
0440         description.constant<double>("HRPPD_COLLECTION_EFFICIENCY");
0441 #endif
0442
0443 #ifdef WITH_IRT2_SUPPORT
0444     // [0,0]: have neither access to G4VSolid nor to G4Material; IRT code does not care; fine;
0445     auto pd = new IRT2::CherenkovPhotonDetector(0, 0);
0446
0447     // FIXME: '0' stands for the unknown (and irrelevant) G4LogicalVolume;
0448     geometry->AddPhotonDetector(cdet, 0, pd);
0449
0450     // Cannot access GEANT shapes in the reconstruction code -> store this value;
0451     pd->SetActiveAreaSize(_HRPPD_ACTIVE_AREA_SIZE_ / mm);
0452
0453     pd->SetGeometricEfficiency(_HRPPD_COLLECTION_EFFICIENCY_);
0454 #endif
0455
0456 #ifdef WITH_IRT2_SUPPORT
0457     {
0458       TVector3 nx(1 * sign, 0, 0), ny(0, -1, 0);
0459
0460       // For now assume it is a unique surface, same for all HRPPDs;
0461       auto surface =
0462           new FlatSurface(sign * (1 / mm) *
0463                               TVector3(0, 0,
0464                                        fvOffset + _FIDUCIAL_VOLUME_LENGTH_ / 2 -
0465                                            _SENSOR_AREA_LENGTH_ + _HRPPD_WINDOW_THICKNESS_ / 2),
0466                           nx, ny);
0467
0468       auto radiator = geometry->AddFlatRadiator(cdet, "QuartzWindow", CherenkovDetector::Downstream,
0469                                                 0, (G4LogicalVolume*)(0x3), 0, surface,
0470                                                 _HRPPD_WINDOW_THICKNESS_ / mm);
0471       radiator->SetAlternativeMaterialName("AirOptical");
0472     }
0473 #endif
0474
0475     // Create a vector of HRPPD XY-coordinates in a separate loop;
0476     std::vector<std::pair<TVector2, bool>> coord;
0477     {
0478       unsigned const hdim              = 9;
0479       const unsigned flags[hdim][hdim] = {
0480           // NB: WYSIWIG fashion; well, it is top/bottom and left/right symmetric;
0481           // clang-format off
0482         {0, 0, 1, 1, 1, 1, 1, 0, 0},
0483         {0, 1, 1, 1, 1, 1, 1, 1, 0},
0484         {1, 1, 1, 1, 1, 1, 1, 1, 1},
0485         {1, 1, 1, 1, 2, 1, 1, 1, 1},
0486         {3, 3, 3, 4, 0, 2, 1, 1, 1},
0487         {1, 1, 1, 1, 2, 1, 1, 1, 1},
0488         {1, 1, 1, 1, 1, 1, 1, 1, 1},
0489         {0, 1, 1, 1, 1, 1, 1, 1, 0},
0490         {0, 0, 1, 1, 1, 1, 1, 0, 0}};
0491       // clang-format on
0492
0493       for (unsigned ix = 0; ix < hdim; ix++) {
0494         double xOffset = (_HRPPD_TILE_SIZE_ + _HRPPD_INSTALLATION_GAP_) * (ix - (hdim - 1) / 2.);
0495
0496         for (unsigned iy = 0; iy < hdim; iy++) {
0497           double yOffset = (_HRPPD_TILE_SIZE_ + _HRPPD_INSTALLATION_GAP_) * (iy - (hdim - 1) / 2.);
0498           unsigned flag  = flags[hdim - iy - 1][ix];
0499
0500           if (!flag)
0501             continue;
0502
0503           double qxOffset = xOffset + (flag >= 3 ? -_HRPPD_CENTRAL_ROW_OFFSET_ : 0.0);
0504           coord.push_back(std::make_pair(TVector2(qxOffset, yOffset), flag % 2));
0505         } //for iy
0506       } //for ix
0507     }
0508
0509     unsigned imod = 0;
0510
0511     //
0512     // It looks like each HRPPD should be a new object rather than a copy of the same volume,
0513     // because otherwise one cannot make photocathodes sensitive (this is done on a PlacedVolume
0514     // level);
0515     //
0516     for (auto xyptr : coord) {
0517       auto& xy = xyptr.first;
0518
0519       uint64_t sensorID = 0x0;
0520
0521       //
0522       // HRPPD container volume
0523       //
0524       Box hrppd_Solid(_HRPPD_TILE_SIZE_ / 2, _HRPPD_TILE_SIZE_ / 2,
0525                       _HRPPD_CONTAINER_VOLUME_HEIGHT_ / 2);
0526       TString hrppdName;
0527       hrppdName.Form("%s-hrppd-%02d", detName.c_str(), imod);
0528       // FIXME: may want to use AirOptical here, but then return a dummy absorber
0529       // layer behind the actual photocathode;
0530       Volume hrppdVol_air(hrppdName.Data(), hrppd_Solid, air);
0531
0532       // A running variable to pack layers one after the other one;
0533       double accu = -_HRPPD_CONTAINER_VOLUME_HEIGHT_ / 2;
0534
0535       //
0536       // Quartz Window
0537       //
0538       Box wnd_Solid(_HRPPD_TILE_SIZE_ / 2, _HRPPD_TILE_SIZE_ / 2, _HRPPD_WINDOW_THICKNESS_ / 2);
0539       TString wndName;
0540       wndName.Form("%s-window-%02d", detName.c_str(), imod);
0541       Volume wndVol(wndName.Data(), wnd_Solid, HRPPD_WindowMat);
0542       wndVol.setVisAttributes(wndVis);
0543       hrppdVol_air.placeVolume(wndVol, Position(0, 0, accu + _HRPPD_WINDOW_THICKNESS_ / 2));
0544
0545       accu += _HRPPD_WINDOW_THICKNESS_;
0546
0547       //
0548       // Photocathode layer (sensitive volume)
0549       //
0550       {
0551         auto pcBox = Box(_HRPPD_ACTIVE_AREA_SIZE_ / 2, _HRPPD_ACTIVE_AREA_SIZE_ / 2,
0552                          _HRPPD_PHOTOCATHODE_THICKNESS_ / 2);
0553         TString pcName;
0554         pcName.Form("%s-photocathode-%02d", detName.c_str(), imod);
0555         Volume pcVol(pcName.Data(), pcBox, HRPPD_pcMat);
0556         pcVol.setSensitiveDetector(sens);
0557
0558         pcVol.setVisAttributes(pcVis);
0559         PlacedVolume pcPV = hrppdVol_air.placeVolume(
0560             pcVol, Position(0.0, 0.0, accu + _HRPPD_PHOTOCATHODE_THICKNESS_ / 2));
0561         {
0562           pcPV.addPhysVolID("hrppd", imod);
0563
0564           // sensor DetElement
0565           sensorID = encodeSensorID(pcPV.volIDs());
0566           TString deName;
0567           deName.Form("%s-sensor-%02d", detName.c_str(), imod);
0568           DetElement pcDE(sdet, deName.Data(), sensorID);
0569           pcDE.setPlacement(pcPV);
0570         }
0571
0572         //
0573         // A fake absorber layer behind the photocathode; FIXME: make sure that reflection
0574         // on the window and photocathode boundary still works correctly (no fake volume as
0575         // in a standalone code);
0576         //
0577         {
0578           TString absName;
0579           absName.Form("%s-absorber-%02d", detName.c_str(), imod);
0580           // Recycle the same pcBox shape; do not mind to use PCB material;
0581           Volume absVol(absName.Data(), pcBox, HRPPD_PCBMat);
0582           hrppdVol_air.placeVolume(absVol, Position(0.0, 0.0,
0583                                                     accu + _HRPPD_PHOTOCATHODE_THICKNESS_ +
0584                                                         _HRPPD_PHOTOCATHODE_THICKNESS_ / 2));
0585         }
0586       }
0587
0588       //
0589       // Ceramic body (sidewall and anode)
0590       //
0591       {
0592         Box cerbox(_HRPPD_TILE_SIZE_ / 2, _HRPPD_TILE_SIZE_ / 2,
0593                    _HRPPD_CERAMIC_BODY_THICKNESS_ / 2);
0594         Box cutbox(_HRPPD_OPEN_AREA_SIZE_ / 2, _HRPPD_OPEN_AREA_SIZE_ / 2,
0595                    _HRPPD_CERAMIC_BODY_THICKNESS_ / 2);
0596
0597         SubtractionSolid ceramic(cerbox, cutbox, Position(0, 0, -_HRPPD_BASEPLATE_THICKNESS_));
0598
0599         TString cerName;
0600         cerName.Form("%s-ceramic-%02d", detName.c_str(), imod);
0601         Volume ceramicVol(cerName.Data(), ceramic, HRPPD_CeramicMat);
0602
0603         ceramicVol.setVisAttributes(bodyVis);
0604         hrppdVol_air.placeVolume(ceramicVol,
0605                                  Position(0.0, 0.0, accu + _HRPPD_CERAMIC_BODY_THICKNESS_ / 2));
0606       }
0607
0608       //
0609       // Effective anode plating layer
0610       //
0611       {
0612         Box plating_solid(_HRPPD_OPEN_AREA_SIZE_ / 2, _HRPPD_OPEN_AREA_SIZE_ / 2,
0613                           _HRPPD_PLATING_LAYER_THICKNESS_ / 2);
0614         TString pltName;
0615         pltName.Form("%s-plating-%02d", detName.c_str(), imod);
0616         Volume platingVol(pltName.Data(), plating_solid, HRPPD_PlatingMat);
0617         // Place somewhere in the middle of the ceramic body gap;
0618         hrppdVol_air.placeVolume(platingVol,
0619                                  Position(0.0, 0.0, accu + _HRPPD_CERAMIC_BODY_THICKNESS_ / 2));
0620       }
0621
0622       //
0623       // Effective MCP layer
0624       //
0625       {
0626         Box mcp_solid(_HRPPD_OPEN_AREA_SIZE_ / 2, _HRPPD_OPEN_AREA_SIZE_ / 2,
0627                       _EFFECTIVE_MCP_THICKNESS_ / 2);
0628         TString mcpName;
0629         mcpName.Form("%s-mcp-%02d", detName.c_str(), imod);
0630         Volume mcpVol(mcpName.Data(), mcp_solid, HRPPD_MCPMat);
0631         hrppdVol_air.placeVolume(mcpVol, Position(0.0, 0.0,
0632                                                   accu + _HRPPD_CERAMIC_BODY_THICKNESS_ / 2 +
0633                                                       _HRPPD_PLATING_LAYER_THICKNESS_ +
0634                                                       _EFFECTIVE_MCP_THICKNESS_ / 2));
0635       }
0636
0637       accu += _HRPPD_CERAMIC_BODY_THICKNESS_;
0638
0639       //
0640       // PCB
0641       //
0642       {
0643         auto _READOUT_PCB_THICKNESS_ = description.constant<double>("READOUT_PCB_THICKNESS");
0644         auto _READOUT_PCB_SIZE_      = description.constant<double>("READOUT_PCB_SIZE");
0645
0646         Box pcb_solid(_READOUT_PCB_SIZE_ / 2, _READOUT_PCB_SIZE_ / 2, _READOUT_PCB_THICKNESS_ / 2);
0647         TString pcbName;
0648         pcbName.Form("%s-pcb-%02d", detName.c_str(), imod);
0649         Volume pcbVol(pcbName.Data(), pcb_solid, HRPPD_PCBMat);
0650         hrppdVol_air.placeVolume(pcbVol, Position(0.0, 0.0, accu + _READOUT_PCB_THICKNESS_ / 2));
0651       }
0652
0653       // Eventually place the whole HRPPD container volume;
0654       double dz =
0655           _FIDUCIAL_VOLUME_LENGTH_ / 2 - _SENSOR_AREA_LENGTH_ + _HRPPD_CONTAINER_VOLUME_HEIGHT_ / 2;
0656       pfRICH_volume.placeVolume(hrppdVol_air, Position(xy.X(), xy.Y(), dz));
0657
0658 #ifdef WITH_IRT2_SUPPORT
0659       {
0660         // Photocathode surface;
0661         double xOffset = xy.X(), yOffset = xy.Y();
0662         auto surface = new FlatSurface(
0663             (1 / mm) *
0664                 TVector3(sign * xOffset, yOffset,
0665                          sign * (fvOffset + _FIDUCIAL_VOLUME_LENGTH_ / 2 - _SENSOR_AREA_LENGTH_ +
0666                                  _HRPPD_WINDOW_THICKNESS_ + _HRPPD_PHOTOCATHODE_THICKNESS_ / 2)),
0667             TVector3(1 * sign, 0, 0), TVector3(0, -1, 0));
0668
0669         {
0670           // '0': pfRICH has no division in sectors (unlike e.g. dRICH);
0671           unsigned sector = 0;
0672
0673           for (unsigned iq = 0; iq < 4; iq++) {
0674             auto irt = pd->AllocateIRT(sector, sensorID);
0675
0676             // Aerogel and acrylic;
0677             if (cdet->m_OpticalBoundaries[CherenkovDetector::Upstream].find(sector) !=
0678                 cdet->m_OpticalBoundaries[CherenkovDetector::Upstream].end())
0679               for (auto boundary : cdet->m_OpticalBoundaries[CherenkovDetector::Upstream][sector])
0680                 irt->AddOpticalBoundary(boundary);
0681
0682             switch (iq) {
0683             case 0:
0684               // Direct hit;
0685               break;
0686             case 1:
0687             case 2:
0688               // Reflection on either inner or outer mirrors;
0689               irt->AddOpticalBoundary(mboundaries[iq - 1]);
0690               break;
0691             case 3:
0692               // Reflection on outer, then on inner mirror; happens at large angles; if the pyramids are
0693               // too high, these photons will undergo more reflections, and cannot be saved;
0694               irt->AddOpticalBoundary(mboundaries[1]);
0695               irt->AddOpticalBoundary(mboundaries[0]);
0696               break;
0697             } //switch
0698
0699             // Fused silica windows;
0700             if (cdet->m_OpticalBoundaries[CherenkovDetector::Downstream].find(sector) !=
0701                 cdet->m_OpticalBoundaries[CherenkovDetector::Downstream].end())
0702               for (auto boundary : cdet->m_OpticalBoundaries[CherenkovDetector::Downstream][sector])
0703                 irt->AddOpticalBoundary(boundary);
0704
0705             // Terminate the optical path;
0706             pd->AddItselfToOpticalBoundaries(irt, surface);
0707           } //for iq
0708         }
0709       }
0710 #endif
0711
0712       imod++;
0713     } //for coord
0714   }
0715
0716   return sdet;
0717 } // createDetector()
0718
0719 // -------------------------------------------------------------------------------------
0720
0721 // clang-format off
0722 DECLARE_DETELEMENT(epic_PFRICH, createDetector)