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 apperture = c.child(_Unicode(apperture));
0092     double app_r        = apperture.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       445.580 * dd4hep::cm;                  //extracted from central_beampipe.xml, line 64
0152   double diameterReduce = 11.0 * dd4hep::cm; //size reduction to avoid overlap with electron pipe
0153   double vacuumDiameterEntrance =
0154       25.792 * dd4hep::cm - diameterReduce; //extracted from central_beampipe.xml, line 64
0155   double vacuumDiameterExit =
0156       17.4 * dd4hep::cm; //15mrad @ entrance to magnet to not overlap electron magnet
0157   double crossingAngle          = -0.025; //radians
0158   double endOfCentralBeamPipe_x = endOfCentralBeamPipe_z * crossingAngle;
0159
0160   //-----------------------------------------------
0161   //calculate gap region center, length, and angle
0162   //-----------------------------------------------
0163
0164   for (int i = 1; i < numMagnets; i++) {
0165
0166     angle_elem_gap.push_back((x_beg[i] - x_end[i - 1]) / (z_beg[i] - z_end[i - 1]));
0167     length_gap.push_back(sqrt(pow(z_beg[i] - z_end[i - 1], 2) + pow(x_beg[i] - x_end[i - 1], 2)));
0168     z_gap.push_back(z_end[i - 1] + 0.5 * length_gap[i - 1] * cos(angle_elem_gap[i - 1]));
0169     x_gap.push_back(x_end[i - 1] + 0.5 * length_gap[i - 1] * sin(angle_elem_gap[i - 1]));
0170   }
0171
0172   //-----------------------------------------------
0173   // fill CutTube storage elements
0174   //-----------------------------------------------
0175
0176   for (int gapIdx = 0; gapIdx < numGaps; gapIdx++) {
0177
0178     inRadius.push_back(0.0);
0179     outRadius.push_back(radii_magnet[gapIdx]);
0180     phi_initial.push_back(0.0);
0181     phi_final.push_back(2 * M_PI);
0182     nxLow.push_back(-(length_gap[gapIdx] / 2.0) *
0183                     sin(rotation_magnet[gapIdx] - angle_elem_gap[gapIdx]));
0184     nyLow.push_back(0.0);
0185     nzLow.push_back(-(length_gap[gapIdx] / 2.0) *
0186                     cos(rotation_magnet[gapIdx] - angle_elem_gap[gapIdx]));
0187     nxHigh.push_back((length_gap[gapIdx] / 2.0) *
0188                      sin(rotation_magnet[gapIdx + 1] - angle_elem_gap[gapIdx]));
0189     nyHigh.push_back(0.0);
0190     nzHigh.push_back((length_gap[gapIdx] / 2.0) *
0191                      cos(rotation_magnet[gapIdx + 1] - angle_elem_gap[gapIdx]));
0192   }
0193
0194   //-----------------------
0195   // inside magnets
0196   //-----------------------
0197
0198   for (int pieceIdx = 0; pieceIdx < numMagnets; pieceIdx++) {
0199
0200     std::string piece_name = Form("MagnetVacuum%d", pieceIdx);
0201
0202     Tube magnetPiece(piece_name, 0.0, radii_magnet[pieceIdx], lengths_magnet[pieceIdx] / 2);
0203     Volume vpiece(piece_name, magnetPiece, m_Vac);
0204     sdet.setAttributes(det, vpiece, x_det.regionStr(), x_det.limitsStr(), vis_name);
0205
0206     assembly.placeVolume(vpiece,
0207                          Transform3D(RotationY(rotation_magnet[pieceIdx]),
0208                                      Position(x_elem_magnet[pieceIdx], y_elem_magnet[pieceIdx],
0209                                               z_elem_magnet[pieceIdx])));
0210   }
0211
0212   //--------------------------
0213   //between magnets
0214   //--------------------------
0215
0216   for (int pieceIdx = numMagnets; pieceIdx < numGaps + numMagnets; pieceIdx++) {
0217
0218     int correctIdx = pieceIdx - numMagnets;
0219
0220     std::string piece_name = Form("GapVacuum%d", correctIdx);
0221
0222     CutTube gapPiece(piece_name, inRadius[correctIdx], outRadius[correctIdx],
0223                      length_gap[correctIdx] / 2, phi_initial[correctIdx], phi_final[correctIdx],
0224                      nxLow[correctIdx], nyLow[correctIdx], nzLow[correctIdx], nxHigh[correctIdx],
0225                      nyHigh[correctIdx], nzHigh[correctIdx]);
0226
0227     Volume vpiece(piece_name, gapPiece, m_Vac);
0228     sdet.setAttributes(det, vpiece, x_det.regionStr(), x_det.limitsStr(), vis_name);
0229
0230     assembly.placeVolume(vpiece, Transform3D(RotationY(angle_elem_gap[correctIdx]),
0231                                              Position(x_gap[correctIdx], 0.0, z_gap[correctIdx])));
0232   }
0233
0234   //--------------------------------------------------------------
0235   //make and place vacuum volume to connect IP beam pipe to B0pf
0236   //--------------------------------------------------------------
0237
0238   if (makeIP_B0pfVacuum) {
0239
0240     double specialGapLength = sqrt(pow(z_beg[0] - endOfCentralBeamPipe_z, 2) +
0241                                    pow(x_beg[0] - endOfCentralBeamPipe_x, 2)) -
0242                               0.1;
0243     double specialGap_z = 0.5 * specialGapLength * cos(crossingAngle) + endOfCentralBeamPipe_z;
0244     double specialGap_x = 0.5 * specialGapLength * sin(crossingAngle) + endOfCentralBeamPipe_x;
0245
0246     std::string piece_name = Form("GapVacuum%d", numGaps + numMagnets);
0247
0248     Cone specialGap(piece_name, specialGapLength / 2, 0.0, vacuumDiameterEntrance / 2, 0.0,
0249                     vacuumDiameterExit / 2);
0250
0251     Volume specialGap_v(piece_name, specialGap, m_Vac);
0252     sdet.setAttributes(det, specialGap_v, x_det.regionStr(), x_det.limitsStr(), vis_name);
0253
0254     assembly.placeVolume(specialGap_v, Transform3D(RotationY(crossingAngle),
0255                                                    Position(specialGap_x, 0.0, specialGap_z)));
0256   }
0257
0258   //----------------------------------------------------
0259
0260   //--------------------------------------------------------------
0261   //make and place vacuum volume after the FF detector array up to end of the world
0262   //--------------------------------------------------------------
0263   if (make_B2pf_EW_Vacuum) {
0264
0265     int pieceIdx           = numMagnets - 1; // last B2PF magnet
0266     std::string piece_name = Form("GapVacuum%d", numGaps + numMagnets + 1);
0267     double endGapLength    = (10000.0 - z_end[pieceIdx]) / cos(rotation_magnet[pieceIdx]);
0268     endGapLength =
0269         endGapLength -
0270         4 * radii_magnet[pieceIdx] *
0271             tan(-rotation_magnet[pieceIdx]); // shift to keep the tube inside the physical volume
0272     double endGap_z = 0.5 * endGapLength * cos(rotation_magnet[pieceIdx]) + z_end[pieceIdx];
0273     double endGap_x = 0.5 * endGapLength * sin(rotation_magnet[pieceIdx]) + x_end[pieceIdx];
0274
0275     Tube vacuum_endWorld(piece_name, 0.0, 4 * radii_magnet[pieceIdx],
0276                          endGapLength / 2); // make larger tube than inner magnet radius
0277     Volume vpiece(piece_name, vacuum_endWorld, m_Vac);
0278     sdet.setAttributes(det, vpiece, x_det.regionStr(), x_det.limitsStr(),
0279                        "InvisibleNoDaughters"); // make invisible instead of AnlBlue
0280
0281     assembly.placeVolume(vpiece, Transform3D(RotationY(rotation_magnet[pieceIdx]),
0282                                              Position(endGap_x, 0.0, endGap_z)));
0283   }
0284   //----------------------------------------------------
0285
0286   pv_assembly = det.pickMotherVolume(sdet).placeVolume(assembly);
0287   pv_assembly.addPhysVolID("system", x_det.id()).addPhysVolID("barrel", 1);
0288   sdet.setPlacement(pv_assembly);
0289   assembly->GetShape()->ComputeBBox();
0290
0291   return sdet;
0292 }
0293
0294 double getRotatedZ(double z, double x, double angle) { return z * cos(angle) - x * sin(angle); }
0295
0296 double getRotatedX(double z, double x, double angle) { return z * sin(angle) + x * cos(angle); }
0297
0298 DECLARE_DETELEMENT(magnetElementInnerVacuum, create_detector)