epic/src/BeamPipeChain_geo.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 - 2025 Dhevan Gangadharan, Simon Gardner
0003
0004 //==========================================================================
0005 //
0006 // Places a chain of beam pipe segments within and between the beamline magnets.
0007 //
0008 // Approximation used for beam pipes in between magnets.
0009 //
0010 //==========================================================================
0011
0012 #include "DD4hep/DetFactoryHelper.h"
0013 #include "DD4hep/Printout.h"
0014 #include "TMath.h"
0015 #include <XML/Helper.h>
0016
0017 using namespace std;
0018 using namespace dd4hep;
0019
0020 static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector /* sens */) {
0021
0022   using namespace ROOT::Math;
0023   xml_det_t x_det = e;
0024   string det_name = x_det.nameStr();
0025   DetElement sdet(det_name, x_det.id());
0026   Assembly assembly(det_name + "_assembly");
0027   Material m_Al     = description.material("Aluminum");
0028   Material m_Vacuum = description.material("Vacuum");
0029   string vis_name   = dd4hep::getAttrOrDefault<std::string>(x_det, _Unicode(vis), "BeamPipeVis");
0030   double thickness  = getAttrOrDefault<double>(x_det, _Unicode(wall_thickness), 0);
0031   double bendRadius = 0.00001 * mm; //Small bend radius to allow the construction of toruses
0032
0033   vector<string> names;
0034   vector<int> ids;
0035   vector<double> xCenters;
0036   vector<double> zCenters;
0037   vector<double> lengths;
0038   vector<double> thetas;
0039   vector<double> rOuters1;
0040   vector<double> rOuters2;
0041   vector<double> bendLengths;
0042
0043   // Grab info for beamline magnets
0044   for (xml_coll_t pipe_coll(x_det, _Unicode(pipe)); pipe_coll; pipe_coll++) { // pipes
0045
0046     xml_comp_t pipe(pipe_coll);
0047
0048     names.push_back(getAttrOrDefault<string>(pipe, _Unicode(name), ""));
0049     ids.push_back(getAttrOrDefault<int>(pipe, _Unicode(id), 0));
0050
0051     // Vectors momentarily filled with zeros for pipes in between magnets
0052     xCenters.push_back(getAttrOrDefault<double>(pipe, _Unicode(xcenter), 0));
0053     zCenters.push_back(getAttrOrDefault<double>(pipe, _Unicode(zcenter), 0));
0054     lengths.push_back(getAttrOrDefault<double>(pipe, _Unicode(length), 0));
0055     thetas.push_back(getAttrOrDefault<double>(pipe, _Unicode(theta), 0));
0056     rOuters1.push_back(getAttrOrDefault<double>(pipe, _Unicode(rout1), 0));
0057     rOuters2.push_back(getAttrOrDefault<double>(pipe, _Unicode(rout2), 0));
0058   }
0059
0060   // Calculate parameters for connecting pipes in between magnets
0061   for (uint pipeN = 0; pipeN < names.size(); pipeN++) {
0062
0063     if (lengths[pipeN] > 0) {
0064       continue;
0065     } // pipe parameters already set to nonzero values
0066     if (pipeN == 0) {
0067       continue;
0068     } // can't create pipe for an empty starting slot
0069     if ((pipeN + 1) == names.size()) {
0070       continue;
0071     } // can't create pipe for an empty end slot
0072
0073     double previousPipeEndX =
0074         xCenters[pipeN - 1] - lengths[pipeN - 1] / 2. * sin(thetas[pipeN - 1]);
0075     double previousPipeEndZ =
0076         zCenters[pipeN - 1] - lengths[pipeN - 1] / 2. * cos(thetas[pipeN - 1]);
0077     double nextPipeStartX = xCenters[pipeN + 1] + lengths[pipeN + 1] / 2. * sin(thetas[pipeN + 1]);
0078     double nextPipeStartZ = zCenters[pipeN + 1] + lengths[pipeN + 1] / 2. * cos(thetas[pipeN + 1]);
0079
0080     double x      = (previousPipeEndX + nextPipeStartX) / 2.;
0081     double z      = (previousPipeEndZ + nextPipeStartZ) / 2.;
0082     double deltaX = previousPipeEndX - nextPipeStartX;
0083     double deltaZ = previousPipeEndZ - nextPipeStartZ;
0084     double l      = sqrt(deltaX * deltaX + deltaZ * deltaZ);
0085     double theta  = atan2(deltaX, deltaZ);
0086
0087     xCenters[pipeN] = x;
0088     zCenters[pipeN] = z;
0089     lengths[pipeN]  = l;
0090     thetas[pipeN]   = theta;
0091     rOuters1[pipeN] = rOuters2[pipeN - 1];
0092     rOuters2[pipeN] = rOuters1[pipeN + 1];
0093   }
0094
0095   // If there is an bend in the pipe, calculate the length reduction of the pipe and joint length
0096   for (uint i = 1; i < thetas.size(); i++) {
0097
0098     // Start at the join between the first two pipes ending at the join between the last two pipes N-1
0099     double bendAngle = thetas[i] - thetas[i - 1];
0100     if (std::abs(bendAngle) < 0.01 * mrad) {
0101       bendLengths.push_back(0);
0102     } else // Correct for tubes, not yet cones so imperfect
0103     {
0104       double bendLength = abs(rOuters1[i] * tan(bendAngle / 2));
0105       bendLengths.push_back(bendLength + bendRadius);
0106     }
0107   }
0108
0109   // Add all pipes to the assembly
0110   for (uint pipeN = 0; pipeN < xCenters.size(); pipeN++) {
0111
0112     double length  = lengths[pipeN];
0113     double theta   = thetas[pipeN];
0114     double xCenter = xCenters[pipeN];
0115     double zCenter = zCenters[pipeN];
0116     double rOuter1 = rOuters1[pipeN];
0117     double rOuter2 = rOuters2[pipeN];
0118
0119     //Change straight pipe length to account for bend
0120     if (pipeN != 0) {
0121       length -= bendLengths[pipeN - 1];
0122       xCenter -= bendLengths[pipeN - 1] / 2 * sin(theta);
0123       zCenter -= bendLengths[pipeN - 1] / 2 * cos(theta);
0124     }
0125     if (pipeN != xCenters.size() - 1) {
0126       length -= bendLengths[pipeN];
0127       xCenter += bendLengths[pipeN] / 2 * sin(theta);
0128       zCenter += bendLengths[pipeN] / 2 * cos(theta);
0129     }
0130
0131     ConeSegment s_tube(length / 2.0, rOuter2 - thickness, rOuter2, rOuter1 - thickness, rOuter1);
0132     ConeSegment s_vacuum(length / 2.0, 0, rOuter2 - thickness, 0, rOuter1 - thickness);
0133
0134     Volume v_tube("v_tube_" + names[pipeN], s_tube, m_Al);
0135     Volume v_vacuum("v_vacuum_" + names[pipeN], s_vacuum, m_Vacuum);
0136
0137     v_tube.setVisAttributes(description.visAttributes(vis_name));
0138
0139     assembly.placeVolume(v_tube, Transform3D(RotationY(theta), Position(xCenter, 0, zCenter)));
0140     auto placed_vacuum = assembly.placeVolume(
0141         v_vacuum, Transform3D(RotationY(theta), Position(xCenter, 0, zCenter)));
0142
0143     DetElement vacuum_element(sdet, names[pipeN] + "_vacuum", ids[pipeN]);
0144     vacuum_element.setPlacement(placed_vacuum);
0145
0146     // // Add joining bend to next pipe
0147     if (pipeN != xCenters.size() - 1 && bendLengths[pipeN] != 0) {
0148
0149       // double bendLength = bendLengths[pipeN];
0150       double bendAngle = thetas[pipeN + 1] - thetas[pipeN];
0151
0152       // Create a torus to join the two pipes
0153       Torus s_bend_solid(rOuter2 + bendRadius, rOuter2 - thickness, rOuter2, 0.0, abs(bendAngle));
0154
0155       //Create a vacuum torus to place inside the solid segment
0156       Torus s_bend_vac(rOuter2 + bendRadius, 0, rOuter2 - thickness, 0, abs(bendAngle));
0157
0158       // //Create the volumes
0159       Volume v_bend_soild("v_bend_solid_" + names[pipeN], s_bend_solid, m_Al);
0160       Volume v_bend_vac("v_bend_vac_" + names[pipeN], s_bend_vac, m_Vacuum);
0161
0162       // Calculate the position and rotation to place bend at the joint
0163       double bendCenterX = xCenter + rOuter2 * cos(theta) - (length / 2) * sin(theta);
0164       double bendCenterZ = zCenter - rOuter2 * sin(theta) - (length / 2) * cos(theta);
0165       double rotation    = pi - theta;
0166       if (bendAngle > 0) {
0167         bendCenterX = xCenter - rOuter2 * cos(theta) - (length / 2) * sin(theta);
0168         bendCenterZ = zCenter + rOuter2 * sin(theta) - (length / 2) * cos(theta);
0169         rotation    = -bendAngle - theta;
0170       }
0171
0172       // Place the bend in the assembly
0173       assembly.placeVolume(v_bend_soild, Transform3D(RotationZYX(rotation, 0, pi / 2),
0174                                                      Position(bendCenterX, 0, bendCenterZ)));
0175       assembly.placeVolume(v_bend_vac, Transform3D(RotationZYX(rotation, 0, pi / 2),
0176                                                    Position(bendCenterX, 0, bendCenterZ)));
0177
0178       // Set vis attributes
0179       v_bend_soild.setVisAttributes(description.visAttributes(vis_name));
0180     }
0181   }
0182
0183   // Final placement
0184   auto pv_assembly =
0185       description.pickMotherVolume(sdet).placeVolume(assembly, Position(0.0, 0.0, 0.0));
0186
0187   sdet.setPlacement(pv_assembly);
0188
0189   assembly->GetShape()->ComputeBBox();
0190
0191   return sdet;
0192 }
0193
0194 DECLARE_DETELEMENT(BeamPipeChain, create_detector)