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