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; //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     // Start at the join between the first two pipes ending at the join between the last two pipes N-1
0098     if (thetas[i - 1] == thetas[i]) {
0099       bendLengths.push_back(0);
0100     } else // Correct for tubes, not yet cones so imperfect
0101     {
0102       double bendAngle  = thetas[i] - thetas[i - 1];
0103       double bendLength = abs(rOuters1[i] * tan(bendAngle / 2));
0104       bendLengths.push_back(bendLength + bendRadius);
0105     }
0106   }
0107
0108   // Add all pipes to the assembly
0109   for (uint pipeN = 0; pipeN < xCenters.size(); pipeN++) {
0110
0111     double length  = lengths[pipeN];
0112     double theta   = thetas[pipeN];
0113     double xCenter = xCenters[pipeN];
0114     double zCenter = zCenters[pipeN];
0115     double rOuter1 = rOuters1[pipeN];
0116     double rOuter2 = rOuters2[pipeN];
0117
0118     //Change straight pipe length to account for bend
0119     if (pipeN != 0) {
0120       length -= bendLengths[pipeN - 1];
0121       xCenter -= bendLengths[pipeN - 1] / 2 * sin(theta);
0122       zCenter -= bendLengths[pipeN - 1] / 2 * cos(theta);
0123     }
0124     if (pipeN != xCenters.size() - 1) {
0125       length -= bendLengths[pipeN];
0126       xCenter += bendLengths[pipeN] / 2 * sin(theta);
0127       zCenter += bendLengths[pipeN] / 2 * cos(theta);
0128     }
0129
0130     ConeSegment s_tube(length / 2.0, rOuter2 - thickness, rOuter2, rOuter1 - thickness, rOuter1);
0131     ConeSegment s_vacuum(length / 2.0, 0, rOuter2 - thickness, 0, rOuter1 - thickness);
0132
0133     Volume v_tube("v_tube_" + names[pipeN], s_tube, m_Al);
0134     Volume v_vacuum("v_vacuum_" + names[pipeN], s_vacuum, m_Vacuum);
0135
0136     v_tube.setVisAttributes(description.visAttributes(vis_name));
0137
0138     assembly.placeVolume(v_tube, Transform3D(RotationY(theta), Position(xCenter, 0, zCenter)));
0139     auto placed_vacuum = assembly.placeVolume(
0140         v_vacuum, Transform3D(RotationY(theta), Position(xCenter, 0, zCenter)));
0141
0142     DetElement vacuum_element(sdet, names[pipeN] + "_vacuum", ids[pipeN]);
0143     vacuum_element.setPlacement(placed_vacuum);
0144
0145     // // Add joining bend to next pipe
0146     if (pipeN != xCenters.size() - 1 && bendLengths[pipeN] != 0) {
0147
0148       // double bendLength = bendLengths[pipeN];
0149       double bendAngle = thetas[pipeN + 1] - thetas[pipeN];
0150
0151       // Create a torus to join the two pipes
0152       Torus s_bend_solid(rOuter2 + bendRadius, rOuter2 - thickness, rOuter2, 0.0, abs(bendAngle));
0153
0154       //Create a vacuum torus to place inside the solid segment
0155       Torus s_bend_vac(rOuter2 + bendRadius, 0, rOuter2 - thickness, 0, abs(bendAngle));
0156
0157       // //Create the volumes
0158       Volume v_bend_soild("v_bend_solid_" + names[pipeN], s_bend_solid, m_Al);
0159       Volume v_bend_vac("v_bend_vac_" + names[pipeN], s_bend_vac, m_Vacuum);
0160
0161       // Calculate the position and rotation to place bend at the joint
0162       double bendCenterX = xCenter + rOuter2 * cos(theta) - (length / 2) * sin(theta);
0163       double bendCenterZ = zCenter - rOuter2 * sin(theta) - (length / 2) * cos(theta);
0164       double rotation    = pi - theta;
0165       if (bendAngle > 0) {
0166         bendCenterX = xCenter - rOuter2 * cos(theta) - (length / 2) * sin(theta);
0167         bendCenterZ = zCenter + rOuter2 * sin(theta) - (length / 2) * cos(theta);
0168         rotation    = -bendAngle - theta;
0169       }
0170
0171       // Place the bend in the assembly
0172       assembly.placeVolume(v_bend_soild, Transform3D(RotationZYX(rotation, 0, pi / 2),
0173                                                      Position(bendCenterX, 0, bendCenterZ)));
0174       assembly.placeVolume(v_bend_vac, Transform3D(RotationZYX(rotation, 0, pi / 2),
0175                                                    Position(bendCenterX, 0, bendCenterZ)));
0176
0177       // Set vis attributes
0178       v_bend_soild.setVisAttributes(description.visAttributes(vis_name));
0179     }
0180   }
0181
0182   // Final placement
0183   auto pv_assembly =
0184       description.pickMotherVolume(sdet).placeVolume(assembly, Position(0.0, 0.0, 0.0));
0185
0186   sdet.setPlacement(pv_assembly);
0187
0188   assembly->GetShape()->ComputeBBox();
0189
0190   return sdet;
0191 }
0192
0193 DECLARE_DETELEMENT(BeamPipeChain, create_detector)