epic/src/magnetVacuumFF.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2023 Alex Jentsch
0003
0004 #include "DD4hep/DetFactoryHelper.h"
0005 #include "DD4hep/Printout.h"
0006 #include <XML/Helper.h>
0007
0008 using namespace std;
0009 using namespace dd4hep;
0010
0011 /** \addtogroup beamline Beamline Instrumentation
0012  */
0013
0014 /** \addtogroup IRChamber Interaction Region Vacuum Chamber.
0015  * \brief Type: **IRChamber**.
0016  * \ingroup beamline
0017  *
0018  *
0019  * \code
0020  *   <detector>
0021  *   </detector>
0022  * \endcode
0023  *
0024  */
0025
0026 static double getRotatedZ(double z, double x, double angle);
0027 static double getRotatedX(double z, double x, double angle);
0028
0029 static Ref_t create_detector(Detector& det, xml_h e, SensitiveDetector /* sens */) {
0030
0031   xml_det_t x_det = e;
0032   string det_name = x_det.nameStr();
0033   DetElement sdet(det_name, x_det.id());
0034   Assembly assembly(det_name + "_assembly");
0035   Material m_Vac  = det.material("Vacuum");
0036   string vis_name = x_det.visStr();
0037
0038   PlacedVolume pv_assembly;
0039
0040   //----------------------------------------------
0041   // Starting point is only the magnet centers,
0042   // lengths, rotations, and radii -->
0043   // everything else calculated internally to
0044   // make it easier to update later.
0045   //----------------------------------------------
0046
0047   bool makeIP_B0pfVacuum   = true; //This is for the special gap location between IP and b0pf
0048   bool make_B2pf_EW_Vacuum = true; //This is for the gap after b2pf
0049
0050   //information for actual FF magnets, with magnet centers as reference
0051   vector<double> radii_magnet;
0052   vector<double> lengths_magnet;
0053   vector<double> rotation_magnet;
0054   vector<double> x_elem_magnet;
0055   vector<double> y_elem_magnet;
0056   vector<double> z_elem_magnet;
0057
0058   //calculated entrance/exit points of FF magnet
0059   vector<double> x_beg;
0060   vector<double> z_beg;
0061   vector<double> x_end;
0062   vector<double> z_end;
0063
0064   //calculated center of gap regions between magnets, rotation, and length
0065   vector<double> angle_elem_gap;
0066   vector<double> z_gap;
0067   vector<double> x_gap;
0068   vector<double> length_gap;
0069
0070   //storage elements for CutTube geometry element used for gaps
0071   vector<double> inRadius;
0072   vector<double> outRadius;
0073   vector<double> nxLow;
0074   vector<double> nyLow;
0075   vector<double> nzLow;
0076   vector<double> nxHigh;
0077   vector<double> nyHigh;
0078   vector<double> nzHigh;
0079   vector<double> phi_initial;
0080   vector<double> phi_final;
0081
0082   for (xml_coll_t c(x_det, _U(element)); c; ++c) {
0083
0084     xml_dim_t pos      = c.child(_U(placement));
0085     double pos_x       = pos.x();
0086     double pos_y       = pos.y();
0087     double pos_z       = pos.z();
0088     double pos_theta   = pos.attr<double>(_U(theta));
0089     xml_dim_t dims     = c.child(_U(dimensions)); //dimensions();
0090     double dim_z       = dims.z();
0091     xml_dim_t aperture = c.child(_Unicode(aperture));
0092     double app_r       = aperture.r();
0093
0094     radii_magnet.push_back(app_r);        // cm
0095     lengths_magnet.push_back(dim_z);      //cm
0096     rotation_magnet.push_back(pos_theta); // radians
0097     x_elem_magnet.push_back(pos_x * dd4hep::cm);
0098     y_elem_magnet.push_back(pos_y * dd4hep::cm);
0099     z_elem_magnet.push_back(pos_z * dd4hep::cm);
0100   }
0101
0102   int numMagnets = radii_magnet.size(); //number of actual FF magnets between IP and FF detectors
0103   int numGaps =
0104       numMagnets -
0105       2; //number of gaps between magnets (excluding the IP to B0pf transition -- special case, and the gao after B1apf)
0106
0107   //-------------------------------------------
0108   // override numbers for the first element -->
0109   // doesn't use the actual B0pf geometry!!!
0110   // -->it's based on the B0 beam pipe
0111   // this needs to be fixed later to read-in
0112   // that beam pipe geometry
0113   //-------------------------------------------
0114
0115   radii_magnet[0]    = 2.9;                 // cm
0116   lengths_magnet[0]  = 120.0;               // cm
0117   rotation_magnet[0] = -0.025;              // radians
0118   x_elem_magnet[0]   = 640.0 * sin(-0.025); // cm
0119   y_elem_magnet[0]   = 0.0;                 // cm
0120   z_elem_magnet[0]   = 640.0 * cos(-0.025); // cm
0121
0122   //-------------------------------------------
0123   //calculate entrance/exit points of magnets
0124   //-------------------------------------------
0125
0126   for (int i = 0; i < numMagnets; i++) {
0127
0128     // need to use the common coordinate system -->
0129     // use x = z, and y = x to make things easier
0130
0131     z_beg.push_back(getRotatedZ(-0.5 * lengths_magnet[i], 0.0, rotation_magnet[i]) +
0132                     z_elem_magnet[i]);
0133     z_end.push_back(getRotatedZ(0.5 * lengths_magnet[i], 0.0, rotation_magnet[i]) +
0134                     z_elem_magnet[i]);
0135     x_beg.push_back(getRotatedX(-0.5 * lengths_magnet[i], 0.0, rotation_magnet[i]) +
0136                     x_elem_magnet[i]);
0137     x_end.push_back(getRotatedX(0.5 * lengths_magnet[i], 0.0, rotation_magnet[i]) +
0138                     x_elem_magnet[i]);
0139   }
0140
0141   //------------------------------------------
0142   // this part is a bit ugly for now -
0143   // it's to make the vacuum volume between the
0144   // end of the IP beam pipe and the beginning of
0145   // beginning of the B0pf magnet
0146   //
0147   // -->the volume will be calculated at the end
0148   //-------------------------------------------
0149
0150   double endOfCentralBeamPipe_z =
0151       494.556 * dd4hep::cm +
0152       2. * dd4hep::mm; //extracted from central_beampipe.xml, line 112 + offset to avoid overlaps
0153   double diameterReduce = 11.0 * dd4hep::cm; //size reduction to avoid overlap with electron pipe
0154   double vacuumDiameterEntrance =
0155       25.792 * dd4hep::cm - diameterReduce; //extracted from central_beampipe.xml, line 64
0156   double vacuumDiameterExit =
0157       17.4 * dd4hep::cm; //15mrad @ entrance to magnet to not overlap electron magnet
0158   double crossingAngle          = -0.025; //radians
0159   double endOfCentralBeamPipe_x = endOfCentralBeamPipe_z * crossingAngle;
0160
0161   //-----------------------------------------------
0162   //calculate gap region center, length, and angle
0163   //-----------------------------------------------
0164
0165   for (int i = 1; i < numMagnets; i++) {
0166
0167     angle_elem_gap.push_back((x_beg[i] - x_end[i - 1]) / (z_beg[i] - z_end[i - 1]));
0168     length_gap.push_back(sqrt(pow(z_beg[i] - z_end[i - 1], 2) + pow(x_beg[i] - x_end[i - 1], 2)));
0169     z_gap.push_back(z_end[i - 1] + 0.5 * length_gap[i - 1] * cos(angle_elem_gap[i - 1]));
0170     x_gap.push_back(x_end[i - 1] + 0.5 * length_gap[i - 1] * sin(angle_elem_gap[i - 1]));
0171   }
0172
0173   //-----------------------------------------------
0174   // fill CutTube storage elements
0175   //-----------------------------------------------
0176
0177   for (int gapIdx = 0; gapIdx < numGaps; gapIdx++) {
0178
0179     inRadius.push_back(0.0);
0180     outRadius.push_back(radii_magnet[gapIdx]);
0181     phi_initial.push_back(0.0);
0182     phi_final.push_back(2 * M_PI);
0183     nxLow.push_back(-(length_gap[gapIdx] / 2.0) *
0184                     sin(rotation_magnet[gapIdx] - angle_elem_gap[gapIdx]));
0185     nyLow.push_back(0.0);
0186     nzLow.push_back(-(length_gap[gapIdx] / 2.0) *
0187                     cos(rotation_magnet[gapIdx] - angle_elem_gap[gapIdx]));
0188     nxHigh.push_back((length_gap[gapIdx] / 2.0) *
0189                      sin(rotation_magnet[gapIdx + 1] - angle_elem_gap[gapIdx]));
0190     nyHigh.push_back(0.0);
0191     nzHigh.push_back((length_gap[gapIdx] / 2.0) *
0192                      cos(rotation_magnet[gapIdx + 1] - angle_elem_gap[gapIdx]));
0193   }
0194
0195   //-----------------------
0196   // inside magnets
0197   //-----------------------
0198
0199   for (int pieceIdx = 0; pieceIdx < numMagnets; pieceIdx++) {
0200
0201     std::string piece_name = Form("MagnetVacuum%d", pieceIdx);
0202
0203     Tube magnetPiece(piece_name, 0.0, radii_magnet[pieceIdx], lengths_magnet[pieceIdx] / 2);
0204     Volume vpiece(piece_name, magnetPiece, m_Vac);
0205     sdet.setAttributes(det, vpiece, x_det.regionStr(), x_det.limitsStr(), vis_name);
0206
0207     assembly.placeVolume(vpiece,
0208                          Transform3D(RotationY(rotation_magnet[pieceIdx]),
0209                                      Position(x_elem_magnet[pieceIdx], y_elem_magnet[pieceIdx],
0210                                               z_elem_magnet[pieceIdx])));
0211   }
0212
0213   //--------------------------
0214   //between magnets
0215   //--------------------------
0216
0217   for (int pieceIdx = numMagnets; pieceIdx < numGaps + numMagnets; pieceIdx++) {
0218
0219     int correctIdx = pieceIdx - numMagnets;
0220
0221     std::string piece_name = Form("GapVacuum%d", correctIdx);
0222
0223     CutTube gapPiece(piece_name, inRadius[correctIdx], outRadius[correctIdx],
0224                      length_gap[correctIdx] / 2, phi_initial[correctIdx], phi_final[correctIdx],
0225                      nxLow[correctIdx], nyLow[correctIdx], nzLow[correctIdx], nxHigh[correctIdx],
0226                      nyHigh[correctIdx], nzHigh[correctIdx]);
0227
0228     Volume vpiece(piece_name, gapPiece, m_Vac);
0229     sdet.setAttributes(det, vpiece, x_det.regionStr(), x_det.limitsStr(), vis_name);
0230
0231     assembly.placeVolume(vpiece, Transform3D(RotationY(angle_elem_gap[correctIdx]),
0232                                              Position(x_gap[correctIdx], 0.0, z_gap[correctIdx])));
0233   }
0234
0235   //--------------------------------------------------------------
0236   //make and place vacuum volume to connect IP beam pipe to B0pf
0237   //--------------------------------------------------------------
0238
0239   if (makeIP_B0pfVacuum) {
0240
0241     double specialGapLength = sqrt(pow(z_beg[0] - endOfCentralBeamPipe_z, 2) +
0242                                    pow(x_beg[0] - endOfCentralBeamPipe_x, 2)) -
0243                               0.1;
0244     double specialGap_z = 0.5 * specialGapLength * cos(crossingAngle) + endOfCentralBeamPipe_z;
0245     double specialGap_x = 0.5 * specialGapLength * sin(crossingAngle) + endOfCentralBeamPipe_x;
0246
0247     std::string piece_name = Form("GapVacuum%d", numGaps + numMagnets);
0248
0249     ConeSegment specialGap(piece_name, specialGapLength / 2, 0.0, vacuumDiameterEntrance / 2, 0.0,
0250                            vacuumDiameterExit / 2, 40 * deg, (360 - 40) * deg);
0251
0252     Volume specialGap_v(piece_name, specialGap, m_Vac);
0253     sdet.setAttributes(det, specialGap_v, x_det.regionStr(), x_det.limitsStr(), vis_name);
0254
0255     assembly.placeVolume(specialGap_v, Transform3D(RotationY(crossingAngle),
0256                                                    Position(specialGap_x, 0.0, specialGap_z)));
0257   }
0258
0259   //----------------------------------------------------
0260
0261   //--------------------------------------------------------------
0262   //make and place vacuum volume after the FF detector array up to end of the world
0263   //--------------------------------------------------------------
0264   if (make_B2pf_EW_Vacuum) {
0265
0266     int pieceIdx           = numMagnets - 1; // last B2PF magnet
0267     std::string piece_name = Form("GapVacuum%d", numGaps + numMagnets + 1);
0268     double endGapLength    = (10000.0 - z_end[pieceIdx]) / cos(rotation_magnet[pieceIdx]);
0269     endGapLength =
0270         endGapLength -
0271         4 * radii_magnet[pieceIdx] *
0272             tan(-rotation_magnet[pieceIdx]); // shift to keep the tube inside the physical volume
0273     double endGap_z = 0.5 * endGapLength * cos(rotation_magnet[pieceIdx]) + z_end[pieceIdx];
0274     double endGap_x = 0.5 * endGapLength * sin(rotation_magnet[pieceIdx]) + x_end[pieceIdx];
0275
0276     Tube vacuum_endWorld(piece_name, 0.0, 4 * radii_magnet[pieceIdx],
0277                          endGapLength / 2); // make larger tube than inner magnet radius
0278     Volume vpiece(piece_name, vacuum_endWorld, m_Vac);
0279     sdet.setAttributes(det, vpiece, x_det.regionStr(), x_det.limitsStr(),
0280                        "InvisibleNoDaughters"); // make invisible instead of AnlBlue
0281
0282     assembly.placeVolume(vpiece, Transform3D(RotationY(rotation_magnet[pieceIdx]),
0283                                              Position(endGap_x, 0.0, endGap_z)));
0284   }
0285   //----------------------------------------------------
0286
0287   pv_assembly = det.pickMotherVolume(sdet).placeVolume(assembly);
0288   pv_assembly.addPhysVolID("system", x_det.id()).addPhysVolID("barrel", 1);
0289   sdet.setPlacement(pv_assembly);
0290   assembly->GetShape()->ComputeBBox();
0291
0292   return sdet;
0293 }
0294
0295 double getRotatedZ(double z, double x, double angle) { return z * cos(angle) - x * sin(angle); }
0296
0297 double getRotatedX(double z, double x, double angle) { return z * sin(angle) + x * cos(angle); }
0298
0299 DECLARE_DETELEMENT(magnetElementInnerVacuum, create_detector)