epic/src/BeamPipeChain_geo.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Dhevan Gangadharan
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 // Right-angled ends with a small air gap to avoid G4 overlap errors
0010 //
0011 //==========================================================================
0012
0013 #include "DD4hep/DetFactoryHelper.h"
0014 #include "DD4hep/Printout.h"
0015 #include "TMath.h"
0016 #include <XML/Helper.h>
0017
0018 using namespace std;
0019 using namespace dd4hep;
0020
0021 static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector /* sens */) {
0022
0023   using namespace ROOT::Math;
0024   xml_det_t x_det = e;
0025   string det_name = x_det.nameStr();
0026   DetElement sdet(det_name, x_det.id());
0027   Assembly assembly(det_name + "_assembly");
0028   Material m_Al     = description.material("Aluminum");
0029   Material m_Vacuum = description.material("Vacuum");
0030   string vis_name   = dd4hep::getAttrOrDefault<std::string>(x_det, _Unicode(vis), "BeamPipeVis");
0031   double thickness  = getAttrOrDefault<double>(x_det, _Unicode(wall_thickness), 0);
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
0042   // Grab info for beamline magnets
0043   for (xml_coll_t pipe_coll(x_det, _Unicode(pipe)); pipe_coll; pipe_coll++) { // pipes
0044
0045     xml_comp_t pipe(pipe_coll);
0046
0047     names.push_back(getAttrOrDefault<string>(pipe, _Unicode(name), ""));
0048     ids.push_back(getAttrOrDefault<int>(pipe, _Unicode(id), 0));
0049
0050     // Vectors momentarily filled with zeros for pipes in between magnets
0051     xCenters.push_back(getAttrOrDefault<double>(pipe, _Unicode(xcenter), 0));
0052     zCenters.push_back(getAttrOrDefault<double>(pipe, _Unicode(zcenter), 0));
0053     lengths.push_back(getAttrOrDefault<double>(pipe, _Unicode(length), 0));
0054     thetas.push_back(getAttrOrDefault<double>(pipe, _Unicode(theta), 0));
0055     rOuters1.push_back(getAttrOrDefault<double>(pipe, _Unicode(rout1), 0));
0056     rOuters2.push_back(getAttrOrDefault<double>(pipe, _Unicode(rout2), 0));
0057   }
0058
0059   // Calculate parameters for connecting pipes in between magnets
0060   for (uint pipeN = 0; pipeN < names.size(); pipeN++) {
0061
0062     if (lengths[pipeN] > 0) {
0063       continue;
0064     } // pipe parameters already set to nonzero values
0065     if (pipeN == 0) {
0066       continue;
0067     } // can't create pipe for an empty starting slot
0068     if ((pipeN + 1) == names.size()) {
0069       continue;
0070     } // can't create pipe for an empty end slot
0071
0072     double x = (xCenters[pipeN - 1] - lengths[pipeN - 1] / 2. * sin(thetas[pipeN - 1]) +
0073                 xCenters[pipeN + 1] + lengths[pipeN + 1] / 2. * sin(thetas[pipeN + 1])) /
0074                2.;
0075     double z = (zCenters[pipeN - 1] - lengths[pipeN - 1] / 2. * cos(thetas[pipeN - 1]) +
0076                 zCenters[pipeN + 1] + lengths[pipeN + 1] / 2. * cos(thetas[pipeN + 1])) /
0077                2.;
0078     double deltaX = (xCenters[pipeN - 1] - lengths[pipeN - 1] / 2. * sin(thetas[pipeN - 1])) -
0079                     (xCenters[pipeN + 1] + lengths[pipeN + 1] / 2. * sin(thetas[pipeN + 1]));
0080     double deltaZ = (zCenters[pipeN - 1] - lengths[pipeN - 1] / 2. * cos(thetas[pipeN - 1])) -
0081                     (zCenters[pipeN + 1] + lengths[pipeN + 1] / 2. * cos(thetas[pipeN + 1]));
0082     double l     = sqrt(pow(deltaX, 2) + pow(deltaZ, 2));
0083     double theta = atan(deltaX / deltaZ);
0084
0085     // Small air gap between connecting and magnet beam pipes to avoid G4 overlap errors
0086     if ((theta != thetas[pipeN - 1]) || (theta != thetas[pipeN + 1])) {
0087       l -= 0.5;
0088     }
0089
0090     xCenters[pipeN] = x;
0091     zCenters[pipeN] = z;
0092     lengths[pipeN]  = l;
0093     thetas[pipeN]   = theta;
0094     rOuters1[pipeN] = rOuters2[pipeN - 1];
0095     rOuters2[pipeN] = rOuters1[pipeN + 1];
0096   }
0097
0098   // Add all pipes to the assembly
0099   for (uint pipeN = 0; pipeN < xCenters.size(); pipeN++) {
0100
0101     ConeSegment s_tube(lengths[pipeN] / 2.0, rOuters2[pipeN] - thickness, rOuters2[pipeN],
0102                        rOuters1[pipeN] - thickness, rOuters1[pipeN]);
0103     ConeSegment s_vacuum(lengths[pipeN] / 2.0, 0, rOuters2[pipeN] - thickness, 0,
0104                          rOuters1[pipeN] - thickness);
0105
0106     Volume v_tube("v_tube_" + names[pipeN], s_tube, m_Al);
0107     Volume v_vacuum("v_vacuum_" + names[pipeN], s_vacuum, m_Vacuum);
0108
0109     v_tube.setVisAttributes(description.visAttributes(vis_name));
0110
0111     assembly.placeVolume(v_tube, Transform3D(RotationY(thetas[pipeN]),
0112                                              Position(xCenters[pipeN], 0, zCenters[pipeN])));
0113     auto placed_vacuum =
0114         assembly.placeVolume(v_vacuum, Transform3D(RotationY(thetas[pipeN]),
0115                                                    Position(xCenters[pipeN], 0, zCenters[pipeN])));
0116
0117     DetElement vacuum_element(sdet, names[pipeN] + "_vacuum", ids[pipeN]);
0118     vacuum_element.setPlacement(placed_vacuum);
0119   }
0120
0121   // Final placement
0122   auto pv_assembly =
0123       description.pickMotherVolume(sdet).placeVolume(assembly, Position(0.0, 0.0, 0.0));
0124
0125   sdet.setPlacement(pv_assembly);
0126
0127   assembly->GetShape()->ComputeBBox();
0128
0129   return sdet;
0130 }
0131
0132 DECLARE_DETELEMENT(BeamPipeChain, create_detector)