File indexing completed on 2025-10-30 07:58:17
0001 
0002 
0003 
0004 
0005 
0006 
0007 #include "DD4hep/DetFactoryHelper.h"
0008 #include "DD4hep/OpticalSurfaces.h"
0009 #include "DD4hep/Printout.h"
0010 #include "DDRec/DetectorData.h"
0011 #include "DDRec/Surface.h"
0012 #include "GeometryHelpers.h"
0013 #include "Math/AxisAngle.h"
0014 #include "Math/Vector3D.h"
0015 #include "Math/VectorUtil.h"
0016 #include "TMath.h"
0017 #include "TString.h"
0018 #include <XML/Helper.h>
0019 
0020 using namespace std;
0021 using namespace dd4hep;
0022 using namespace dd4hep::rec;
0023 
0024 using Placements = vector<PlacedVolume>;
0025 
0026 static Ref_t createDetector(Detector& description, xml::Handle_t e, SensitiveDetector sens) {
0027   xml_det_t x_det = e;
0028   Material air    = description.material("AirOptical");
0029   string det_name = x_det.nameStr();
0030   DetElement sdet(det_name, x_det.id());
0031   Assembly assembly(det_name);
0032   sens.setType("tracker");
0033   OpticalSurfaceManager surfMgr = description.surfaceManager();
0034 
0035   bool projective = getAttrOrDefault(x_det, _Unicode(projective), false);
0036   bool reflect    = x_det.reflect(true);
0037 
0038   PlacedVolume pv;
0039 
0040   map<string, Volume> modules;
0041   map<string, Placements> sensitives;
0042   map<string, Volume> module_assemblies;
0043   std::map<std::string, DetElement> module_assembly_delements;
0044 
0045   int n_sensor = 1;
0046 
0047   
0048   xml::Component dims = x_det.dimensions();
0049   auto rmin           = dims.rmin();
0050   auto rmax           = dims.rmax();
0051   auto length         = dims.length();
0052   auto zmin           = dims.zmin();
0053   auto zpos           = zmin + length / 2;
0054 
0055   
0056   Tube envShape(rmin, rmax, length / 2., 0., 2 * M_PI);
0057   Volume envVol("MRICH_Envelope", envShape, air);
0058   envVol.setVisAttributes(description.visAttributes(x_det.visStr()));
0059   if (x_det.hasChild(_Unicode(envelope))) {
0060     xml_comp_t x_envelope = x_det.child(_Unicode(envelope));
0061     double thickness      = x_envelope.thickness();
0062     Material material     = description.material(x_envelope.materialStr());
0063     Tube envInsideShape(rmin + thickness, rmax - thickness, length / 2. - thickness);
0064     SubtractionSolid envShellShape(envShape, envInsideShape);
0065     Volume envShell("MRICH_Envelope_Inside", envShellShape, material);
0066     envVol.placeVolume(envShell);
0067   }
0068 
0069   
0070   xml_comp_t x_mod  = x_det.child(_U(module));
0071   string mod_name   = x_mod.nameStr();
0072   double mod_width  = getAttrOrDefault(x_mod, _U(width), 130.0 * mm);
0073   double mod_height = getAttrOrDefault(x_mod, _U(height), 130.0 * mm);
0074   double mod_length = getAttrOrDefault(x_mod, _U(length), 130.0 * mm);
0075 
0076   
0077   Box m_solid(mod_width / 2.0, mod_height / 2.0, mod_length / 2.0);
0078   Volume m_volume(mod_name, m_solid, air);
0079   m_volume.setVisAttributes(description.visAttributes(x_mod.visStr()));
0080   DetElement mod_de(mod_name + std::string("_mod_") + std::to_string(1), 1);
0081   double z_placement = -mod_length / 2.0;
0082 
0083   
0084   if (x_mod.hasChild(_Unicode(frame))) {
0085     xml_comp_t x_frame     = x_mod.child(_Unicode(frame));
0086     double frame_thickness = getAttrOrDefault(x_frame, _U(thickness), 2.0 * mm);
0087     Box frame_inside(mod_width / 2.0 - frame_thickness, mod_height / 2.0 - frame_thickness,
0088                      mod_length / 2.0 - frame_thickness);
0089     SubtractionSolid frame_solid(m_solid, frame_inside);
0090     Material frame_mat = description.material(x_frame.materialStr());
0091     Volume frame_vol(mod_name + "_frame", frame_solid, frame_mat);
0092     auto frame_vis = getAttrOrDefault<std::string>(x_frame, _U(vis), std::string("GrayVis"));
0093     frame_vol.setVisAttributes(description.visAttributes(frame_vis));
0094     
0095     z_placement += frame_thickness / 2.0;
0096     
0097     m_volume.placeVolume(frame_vol);
0098     
0099     z_placement += frame_thickness / 2.0;
0100   }
0101 
0102   
0103   if (x_mod.hasChild(_Unicode(aerogel))) {
0104     xml_comp_t x_aerogel  = x_mod.child(_Unicode(aerogel));
0105     double aerogel_width  = getAttrOrDefault(x_aerogel, _U(width), 130.0 * mm);
0106     double aerogel_length = getAttrOrDefault(x_aerogel, _U(length), 130.0 * mm);
0107     Material aerogel_mat  = description.material(x_aerogel.materialStr());
0108     auto aerogel_vis =
0109         getAttrOrDefault<std::string>(x_aerogel, _U(vis), std::string("InvisibleWithDaughters"));
0110 
0111     xml_comp_t x_aerogel_frame = x_aerogel.child(_Unicode(frame));
0112     double foam_thickness      = getAttrOrDefault(x_aerogel_frame, _U(thickness), 2.0 * mm);
0113     Material foam_mat          = description.material(x_aerogel_frame.materialStr());
0114     auto foam_vis = getAttrOrDefault<std::string>(x_aerogel_frame, _U(vis), std::string("RedVis"));
0115 
0116     
0117     Box foam_box(aerogel_width / 2.0 + foam_thickness, aerogel_width / 2.0 + foam_thickness,
0118                  (aerogel_length + foam_thickness) / 2.0);
0119     Box foam_sub_box(aerogel_width / 2.0, aerogel_width / 2.0,
0120                      (aerogel_length + foam_thickness) / 2.0);
0121     SubtractionSolid foam_frame_solid(foam_box, foam_sub_box, Position(0, 0, foam_thickness));
0122     Volume foam_vol(mod_name + "_aerogel_frame", foam_frame_solid, foam_mat);
0123     foam_vol.setVisAttributes(description.visAttributes(foam_vis));
0124 
0125     
0126     Box aerogel_box(aerogel_width / 2.0, aerogel_width / 2.0, (aerogel_length) / 2.0);
0127     Volume aerogel_vol(mod_name + "_aerogel", aerogel_box, aerogel_mat);
0128     aerogel_vol.setVisAttributes(description.visAttributes(aerogel_vis));
0129 
0130     
0131     z_placement += (aerogel_length + foam_thickness) / 2.0;
0132     
0133     pv = m_volume.placeVolume(foam_vol, Position(0, 0, z_placement));
0134     
0135     z_placement += foam_thickness / 2.0;
0136     pv = m_volume.placeVolume(aerogel_vol, Position(0, 0, z_placement));
0137     DetElement aerogel_de(mod_de, mod_name + std::string("_aerogel_de") + std::to_string(1), 1);
0138     aerogel_de.setPlacement(pv);
0139     
0140     z_placement += aerogel_length / 2.0;
0141 
0142     
0143     auto aerogel_surf = surfMgr.opticalSurface(dd4hep::getAttrOrDefault<std::string>(
0144         x_aerogel, _Unicode(surface), "MRICH_AerogelOpticalSurface"));
0145     SkinSurface skin_surf(description, aerogel_de, Form("MRICH_aerogel_skin_surface_%d", 1),
0146                           aerogel_surf, aerogel_vol);
0147     skin_surf.isValid();
0148   }
0149 
0150   
0151   if (x_mod.hasChild(_Unicode(lens))) {
0152     xml_comp_t x_lens = x_mod.child(_Unicode(lens));
0153 
0154     
0155     
0156     
0157     auto lens_vis       = getAttrOrDefault<std::string>(x_lens, _U(vis), std::string("AnlBlue"));
0158     double groove_pitch = getAttrOrDefault(x_lens, _Unicode(pitch), 0.2 * mm); 
0159     double lens_f       = getAttrOrDefault(x_lens, _Unicode(focal_length), 6.0 * 2.54 * cm);
0160     double eff_diameter = getAttrOrDefault(x_lens, _Unicode(effective_diameter), 152.4 * mm);
0161     double lens_width   = getAttrOrDefault(x_lens, _Unicode(width), 6.7 * 2.54 * cm);
0162     double center_thickness =
0163         getAttrOrDefault(x_lens, _U(thickness), 0.068 * 2.54 * cm); 
0164 
0165     double n_acrylic      = 1.49;
0166     double lens_curvature = 1.0 / (lens_f * (n_acrylic - 1.0)); 
0167     double full_ring_rmax = std::min(eff_diameter / 2.0, lens_width / 2.0);
0168 
0169     double N_grooves        = std::ceil((full_ring_rmax) / groove_pitch);
0170     double groove_last_rmin = (N_grooves - 1) * groove_pitch;
0171     double groove_last_rmax = N_grooves * groove_pitch;
0172 
0173     auto groove_sagitta = [&](double r) { return lens_curvature * std::pow(r, 2) / (1.0 + 1.0); };
0174     double lens_thickness =
0175         groove_sagitta(groove_last_rmax) - groove_sagitta(groove_last_rmin) + center_thickness;
0176 
0177     Material lens_mat = description.material(x_lens.materialStr());
0178     Box lens_box(lens_width / 2.0, lens_width / 2.0, (center_thickness) / 2.0);
0179     SubtractionSolid flat_lens(lens_box, Tube(0.0, full_ring_rmax, 2 * center_thickness));
0180 
0181     Assembly lens_vol(mod_name + "_lens");
0182     Volume flatpart_lens_vol("flatpart_lens", flat_lens, lens_mat);
0183     lens_vol.placeVolume(flatpart_lens_vol);
0184 
0185     int i_groove       = 0;
0186     double groove_rmax = groove_pitch;
0187     double groove_rmin = 0;
0188 
0189     while (groove_rmax <= full_ring_rmax) {
0190       double dZ = groove_sagitta(groove_rmax) - groove_sagitta(groove_rmin);
0191       Polycone groove_solid(
0192           0, 2.0 * M_PI, {groove_rmin, groove_rmin, groove_rmin},
0193           {groove_rmax, groove_rmax, groove_rmin},
0194           {-lens_thickness / 2.0, lens_thickness / 2.0 - dZ, lens_thickness / 2.0});
0195       Volume lens_groove_vol("lens_groove_" + std::to_string(i_groove), groove_solid, lens_mat);
0196       lens_vol.placeVolume(lens_groove_vol);
0197 
0198       i_groove++;
0199       groove_rmin = (i_groove)*groove_pitch;
0200       groove_rmax = (i_groove + 1) * groove_pitch;
0201     }
0202 
0203     lens_vol.setVisAttributes(description.visAttributes(lens_vis));
0204 
0205     
0206     z_placement += lens_thickness / 2.0;
0207     
0208     pv = m_volume.placeVolume(lens_vol, Position(0, 0, z_placement));
0209     DetElement lens_de(mod_de, mod_name + std::string("_lens_de") + std::to_string(1), 1);
0210     lens_de.setPlacement(pv);
0211     
0212     z_placement += lens_thickness / 2.0;
0213 
0214     
0215     auto lens_surf = surfMgr.opticalSurface(dd4hep::getAttrOrDefault<std::string>(
0216         x_lens, _Unicode(surface), "MRICH_LensOpticalSurface"));
0217     SkinSurface skin_surf(description, lens_de, Form("MRichFresnelLens_skin_surface_%d", 1),
0218                           lens_surf, lens_vol);
0219     skin_surf.isValid();
0220   }
0221 
0222   
0223   if (x_mod.hasChild(_Unicode(space))) {
0224     xml_comp_t x_space = x_mod.child(_Unicode(space));
0225     z_placement += getAttrOrDefault(x_space, _U(thickness), 0.0 * mm);
0226   }
0227 
0228   
0229   if (x_mod.hasChild(_Unicode(mirror))) {
0230     xml_comp_t x_mirror  = x_mod.child(_Unicode(mirror));
0231     auto mirror_vis      = getAttrOrDefault<std::string>(x_mirror, _U(vis), std::string("AnlGray"));
0232     double mirror_x1     = getAttrOrDefault(x_mirror, _U(x1), 100.0 * mm);
0233     double mirror_x2     = getAttrOrDefault(x_mirror, _U(x2), 80.0 * mm);
0234     double mirror_length = getAttrOrDefault(x_mirror, _U(length), 130.0 * mm);
0235     double mirror_thickness = getAttrOrDefault(x_mirror, _U(thickness), 2.0 * mm);
0236     double outer_x1         = (mirror_x1 + mirror_thickness) / 2.0;
0237     double outer_x2         = (mirror_x2 + mirror_thickness) / 2.0;
0238     Trd2 outer_mirror_trd(outer_x1, outer_x2, outer_x1, outer_x2, mirror_length / 2.0);
0239     Trd2 inner_mirror_trd(mirror_x1 / 2.0, mirror_x2 / 2.0, mirror_x1 / 2.0, mirror_x2 / 2.0,
0240                           mirror_length / 2.0 + 0.1 * mm);
0241     SubtractionSolid mirror_solid(outer_mirror_trd, inner_mirror_trd);
0242     Material mirror_mat = description.material(x_mirror.materialStr());
0243     Volume mirror_vol(mod_name + "_mirror", mirror_solid, mirror_mat);
0244 
0245     
0246     z_placement += mirror_length / 2.0;
0247     
0248     pv = m_volume.placeVolume(mirror_vol, Position(0, 0, z_placement));
0249     DetElement mirror_de(mod_de, mod_name + std::string("_mirror_de") + std::to_string(1), 1);
0250     mirror_de.setPlacement(pv);
0251     
0252     z_placement += mirror_length / 2.0;
0253 
0254     
0255     auto mirror_surf = surfMgr.opticalSurface(dd4hep::getAttrOrDefault<std::string>(
0256         x_mirror, _Unicode(surface), "MRICH_MirrorOpticalSurface"));
0257     SkinSurface skin_surf(description, mirror_de, Form("MRICH_mirror_skin_surface_%d", 1),
0258                           mirror_surf, mirror_vol);
0259     skin_surf.isValid();
0260   }
0261 
0262   
0263   if (x_mod.hasChild(_Unicode(photodet))) {
0264     xml_comp_t x_photodet = x_mod.child(_Unicode(photodet));
0265     auto photodet_vis = getAttrOrDefault<std::string>(x_photodet, _U(vis), std::string("AnlRed"));
0266     double photodet_width     = getAttrOrDefault(x_photodet, _U(width), 130.0 * mm);
0267     double photodet_thickness = getAttrOrDefault(x_photodet, _U(thickness), 2.0 * mm);
0268     Material photodet_mat     = description.material(x_photodet.materialStr());
0269     Box window_box(photodet_width / 2.0, photodet_width / 2.0, photodet_thickness / 2.0);
0270     Volume window_vol(mod_name + "_window", window_box, photodet_mat);
0271 
0272     
0273     z_placement += photodet_thickness / 2.0;
0274     
0275     pv = m_volume.placeVolume(window_vol, Position(0, 0, z_placement));
0276     DetElement comp_de(mod_de, mod_name + std::string("_sensor_de_") + std::to_string(1), 1);
0277     comp_de.setPlacement(pv);
0278     
0279     z_placement += photodet_thickness / 2.0;
0280 
0281     
0282     pv.addPhysVolID("sensor", n_sensor);
0283     window_vol.setSensitiveDetector(sens);
0284     sensitives[mod_name].push_back(pv);
0285     ++n_sensor;
0286 
0287     
0288     
0289     
0290     
0291     
0292     
0293     
0294 
0295     
0296     int i_layer = 1;
0297     for (xml_coll_t li(x_photodet, _Unicode(layer)); li; ++li) {
0298       xml_comp_t x_layer     = li;
0299       Material layer_mat     = description.material(x_layer.materialStr());
0300       double layer_thickness = x_layer.thickness();
0301       Box layer_box(photodet_width / 2.0, photodet_width / 2.0, layer_thickness / 2.0);
0302       Volume layer_vol(mod_name + "_layer_" + std::to_string(i_layer), layer_box, layer_mat);
0303 
0304       
0305       z_placement += layer_thickness / 2.0;
0306       
0307       pv = m_volume.placeVolume(layer_vol, Position(0, 0, z_placement));
0308       DetElement layer_de(mod_de, mod_name + std::string("_layer_de_") + std::to_string(i_layer),
0309                           1);
0310       layer_de.setPlacement(pv);
0311       
0312       z_placement += layer_thickness / 2.0;
0313 
0314       i_layer++;
0315     }
0316   }
0317 
0318   
0319   
0320   
0321   
0322   
0323   
0324   
0325   
0326   
0327   
0328   
0329   
0330   
0331 
0332   modules[mod_name]                   = m_volume;
0333   module_assembly_delements[mod_name] = mod_de;
0334   
0335 
0336   
0337   auto points = epic::geo::fillSquares({0., 0.}, mod_width, rmin, rmax);
0338 
0339   
0340   auto mod_v = modules[mod_name];
0341 
0342   
0343   std::vector<std::tuple<double, double, double>> positions;
0344   for (xml_coll_t x_positions_i(x_det, _Unicode(positions)); x_positions_i; ++x_positions_i) {
0345     xml_comp_t x_positions = x_positions_i;
0346     for (xml_coll_t x_position_i(x_positions, _U(position)); x_position_i; ++x_position_i) {
0347       xml_comp_t x_position = x_position_i;
0348       positions.push_back(std::make_tuple(x_positions.scale() * x_position.x() * mm,
0349                                           x_positions.scale() * x_position.y() * mm,
0350                                           -x_positions.z0()));
0351     }
0352   }
0353   
0354   if (positions.empty()) {
0355     for (double x = mod_width / 2.0; x < rmax - mod_width / 2.0; x += mod_width) {
0356       for (double y = mod_width / 2.0; y < rmax - mod_width / 2.0; y += mod_width) {
0357         if (pow(x + mod_width / 2.0, 2) + pow(y + mod_width / 2.0, 2) > rmax * rmax)
0358           continue;
0359         if (pow(x - mod_width / 2.0, 2) + pow(y - mod_width / 2.0, 2) < rmin * rmin)
0360           continue;
0361         positions.push_back(std::make_tuple(x, y, 0));
0362       }
0363     }
0364   }
0365 
0366   
0367   int i_mod = 1; 
0368   for (auto& position : positions) {
0369 
0370     
0371     double position_x  = std::get<0>(position);
0372     double position_y  = std::get<1>(position);
0373     double position_z0 = std::get<2>(position);
0374 
0375     
0376     for (auto& p : decltype(positions){{position_x, position_y, position_z0},
0377                                        {position_y, -position_x, position_z0},
0378                                        {-position_x, -position_y, position_z0},
0379                                        {-position_y, position_x, position_z0}}) {
0380 
0381       
0382       double x  = std::get<0>(p);
0383       double y  = std::get<1>(p);
0384       double z0 = std::get<2>(p);
0385 
0386       Transform3D tr;
0387       if (projective) {
0388         double rotAngX = atan(y / z0);
0389         double rotAngY = -1. * atan(x / z0);
0390         tr             = Translation3D(x, y, 0) * RotationX(rotAngX) * RotationY(rotAngY);
0391       } else {
0392         tr = Translation3D(x, y, 0) * RotationX(0);
0393       }
0394 
0395       
0396       pv = envVol.placeVolume(mod_v, tr);
0397       pv.addPhysVolID("module", i_mod);
0398 
0399       auto mod_det_element =
0400           module_assembly_delements[mod_name].clone(mod_name + "__" + std::to_string(i_mod));
0401       mod_det_element.setPlacement(pv);
0402       sdet.add(mod_det_element);
0403 
0404       i_mod++;
0405     }
0406   }
0407 
0408   
0409   if (x_det.hasChild(_Unicode(layer))) {
0410     xml_comp_t x_layer     = x_det.child(_Unicode(layer));
0411     double layer_thickness = x_layer.thickness();
0412     Material layer_mat     = description.material(x_layer.materialStr());
0413     Tube frameShape(rmin, rmax, layer_thickness / 2., 0., 2 * M_PI);
0414     Volume frameVol("MRICH_Frame", frameShape, layer_mat);
0415     pv = envVol.placeVolume(frameVol, Position(0, 0, (length - layer_thickness) / 2.0));
0416   }
0417 
0418   
0419   Volume motherVol = description.pickMotherVolume(sdet);
0420   if (reflect) {
0421     pv = motherVol.placeVolume(envVol, Transform3D(RotationZYX(0, M_PI, 0), Position(0, 0, -zpos)));
0422   } else {
0423     pv = motherVol.placeVolume(envVol, Transform3D(RotationZYX(0, 0, 0), Position(0, 0, +zpos)));
0424   }
0425   pv.addPhysVolID("system", x_det.id());
0426   sdet.setPlacement(pv);
0427   return sdet;
0428 }
0429 
0430 
0431 
0432 
0433 
0434 
0435 
0436 
0437 
0438 
0439 
0440 
0441 
0442 
0443 
0444 
0445 
0446 
0447 
0448 
0449 
0450 
0451 
0452 
0453 
0454 
0455 
0456 
0457 
0458 
0459 
0460 
0461 
0462 
0463 
0464 
0465 
0466 
0467 
0468 
0469 
0470 
0471 
0472 
0473 
0474 
0475 DECLARE_DETELEMENT(epic_MRICH, createDetector)