macros/common/G4_BeamLine.C

0001 #ifndef MACRO_G4BEAMLINE_C
0002 #define MACRO_G4BEAMLINE_C
0003
0004 #include <GlobalVariables.C>
0005
0006 #include <G4_Pipe.C>
0007
0008 #include <g4detectors/BeamLineMagnetSubsystem.h>
0009 #include <g4detectors/PHG4BlockSubsystem.h>
0010 #include <g4detectors/PHG4ConeSubsystem.h>
0011 #include <g4detectors/PHG4CylinderSubsystem.h>
0012 #include <g4detectors/PHG4DetectorSubsystem.h>
0013 #include <g4detectors/PHG4ZDCSubsystem.h>
0014
0015 #include <g4main/PHG4Reco.h>
0016
0017 #include <Rtypes.h>  // resolves R__LOAD_LIBRARY for clang-tidy
0018
0019 #include <algorithm>
0020 #include <fstream>
0021 #include <set>
0022
0023 R__LOAD_LIBRARY(libg4detectors.so)
0024
0025 float PosFlip(float pos);
0026 float AngleFlip(float angle);
0027 float MagFieldFlip(float Bfield);
0028
0029 // This creates the Enable Flag to be used in the main steering macro
0030 namespace Enable
0031 {
0032   bool BEAMLINE = false;
0033   bool BEAMLINE_ABSORBER = false;
0034   bool BEAMLINE_BLACKHOLE = false;
0035   bool BEAMLINE_OVERLAPCHECK = false;
0036   int BEAMLINE_VERBOSITY = 0;
0037
0038 }  // namespace Enable
0039
0040 namespace G4BEAMLINE
0041 {
0042   // the beampipes seem to add 2 no_overlaps - needs to be looked at
0043   // but this z position takes care of our current overlap issues
0044   double starting_z{0};
0045   double enclosure_z_max{0};
0046   double enclosure_r_max = 30.;  // 30cm radius to cover magnets
0047   double enclosure_center{0};
0048   double skin_thickness = 0.;  // if center of magnet iron is black hole - thickness of Fe surrounding it
0049   int pipe_id_offset = 100;
0050   int roman_pot_pipe_id_offset = 200;
0051   PHG4CylinderSubsystem *ForwardBeamLineEnclosure(nullptr);
0052   PHG4CylinderSubsystem *BackwardBeamLineEnclosure(nullptr);
0053
0054 }  // namespace G4BEAMLINE
0055
0056 void BeamLineInit()
0057 {
0058   int verbosity = std::max(Enable::VERBOSITY, Enable::BEAMLINE_VERBOSITY);
0059
0060   // this is a cheap trick to deal with a shiftd beam pipe, we have now gaps between the end of the sPHENIX beam pipe
0061   // and the forward/backward vacuum enclosures. Sadly the magnet and ZDC positions are relative to this volume
0062   // it's not clear if we have two different sizes forward/backward how this would affect the positioning
0063   // of the beam line magnets and the zdcs. This here should preserve them.
0064   G4BEAMLINE::starting_z = std::max(std::abs(G4PIPE::max_z_north), std::abs(G4PIPE::max_z_south));
0065   if (Enable::MVTX_APPLYMISALIGNMENT && !Enable::PIPE)
0066   {
0067     if (verbosity > 0)
0068     {
0069       std::cout << "G4BEAMLINE: pipe misaligned but pipe disabled - applying zshift" << std::endl;
0070     }
0071     G4BEAMLINE::starting_z += G4PIPE::pipe_zshift_2024;
0072   }
0073   G4BEAMLINE::enclosure_z_max = 2050. + (800 - G4BEAMLINE::starting_z);
0074   G4BEAMLINE::enclosure_center = 0.5 * (G4BEAMLINE::starting_z + G4BEAMLINE::enclosure_z_max);
0075
0076   BlackHoleGeometry::min_z = std::min(BlackHoleGeometry::min_z, -G4BEAMLINE::enclosure_z_max);
0077   BlackHoleGeometry::max_z = std::max(BlackHoleGeometry::max_z, G4BEAMLINE::enclosure_z_max);
0078   BlackHoleGeometry::max_radius = std::max(BlackHoleGeometry::max_radius, G4BEAMLINE::enclosure_r_max);
0079 }
0080
0081 void BeamLineDefineMagnets(PHG4Reco *g4Reco)
0082 {
0083   bool overlapCheck = Enable::OVERLAPCHECK || Enable::BEAMLINE_OVERLAPCHECK;
0084   bool AbsorberActive = Enable::ABSORBER || Enable::BEAMLINE_ABSORBER;
0085   int verbosity = std::max(Enable::VERBOSITY, Enable::BEAMLINE_VERBOSITY);
0086
0087   G4BEAMLINE::ForwardBeamLineEnclosure = new PHG4CylinderSubsystem("ForwardBeamLineEnclosure");
0088   G4BEAMLINE::ForwardBeamLineEnclosure->set_double_param("place_z", G4BEAMLINE::enclosure_center);
0089   G4BEAMLINE::ForwardBeamLineEnclosure->set_double_param("radius", 0);
0090   G4BEAMLINE::ForwardBeamLineEnclosure->set_double_param("thickness", G4BEAMLINE::enclosure_r_max);
0091   G4BEAMLINE::ForwardBeamLineEnclosure->set_double_param("length", G4BEAMLINE::enclosure_z_max - G4BEAMLINE::starting_z);
0092   G4BEAMLINE::ForwardBeamLineEnclosure->set_string_param("material", "G4_Galactic");
0093   G4BEAMLINE::ForwardBeamLineEnclosure->set_color(.5, .5, .5, 0.2);
0094   G4BEAMLINE::ForwardBeamLineEnclosure->OverlapCheck(overlapCheck);
0095   if (verbosity)
0096   {
0097     G4BEAMLINE::ForwardBeamLineEnclosure->Verbosity(verbosity);
0098   }
0099   g4Reco->registerSubsystem(G4BEAMLINE::ForwardBeamLineEnclosure);
0100
0101   G4BEAMLINE::BackwardBeamLineEnclosure = new PHG4CylinderSubsystem("BackwardBeamLineEnclosure");
0102   G4BEAMLINE::BackwardBeamLineEnclosure->set_double_param("place_z", -G4BEAMLINE::enclosure_center);
0103   G4BEAMLINE::BackwardBeamLineEnclosure->set_double_param("radius", 0);
0104   G4BEAMLINE::BackwardBeamLineEnclosure->set_double_param("thickness", G4BEAMLINE::enclosure_r_max);  // This is intentionally made large 25cm radius
0105   G4BEAMLINE::BackwardBeamLineEnclosure->set_double_param("length", G4BEAMLINE::enclosure_z_max - G4BEAMLINE::starting_z);
0106   G4BEAMLINE::BackwardBeamLineEnclosure->set_string_param("material", "G4_Galactic");
0107   G4BEAMLINE::BackwardBeamLineEnclosure->set_color(.5, .5, .5, 0.2);
0108   G4BEAMLINE::BackwardBeamLineEnclosure->OverlapCheck(overlapCheck);
0109   if (verbosity)
0110   {
0111     G4BEAMLINE::BackwardBeamLineEnclosure->Verbosity(verbosity);
0112   }
0113   g4Reco->registerSubsystem(G4BEAMLINE::BackwardBeamLineEnclosure);
0114
0115   std::string magFile = std::string(getenv("CALIBRATIONROOT")) + "/Beam/D0DXMagnets.dat";
0116
0117   // if you insert numbers it only displays those magnets, do not comment out the set declaration
0118   std::set<int> magnetlist;
0119   // magnetlist.insert(7);
0120
0121   BeamLineMagnetSubsystem *bl = nullptr;
0122   std::ifstream infile(magFile);
0123   if (infile.is_open())
0124   {
0125     double biggest_z = 0.;
0126     int imagnet = 0;
0127     std::string line;
0128     while (std::getline(infile, line))
0129     {
0130       if (!line.compare(0, 1, "B") ||
0131           !line.compare(0, 1, "Q") ||
0132           !line.compare(0, 1, "S"))
0133       {
0134         std::istringstream iss(line);
0135         std::string magname;
0136         double x;
0137         double y;
0138         double z;
0139         double inner_radius_zin;
0140         double inner_radius_zout;
0141         double outer_magnet_diameter;
0142         double length;
0143         double angle;
0144         double dipole_field_x;
0145         double fieldgradient;
0146         if (!(iss >> magname >> x >> y >> z >> inner_radius_zin >> inner_radius_zout >> outer_magnet_diameter >> length >> angle >> dipole_field_x >> fieldgradient))
0147         {
0148           std::cout << "could not decode " << line << std::endl;
0149           gSystem->Exit(1);
0150         }
0151         else
0152         {
0153           //------------------------
0154
0155           std::string magtype;
0156           if (inner_radius_zin != inner_radius_zout)
0157           {
0158             std::cout << "inner radius at front of magnet " << inner_radius_zin
0159                       << " not equal radius at back of magnet " << inner_radius_zout
0160                       << " needs change in code (replace tube by cone for beamline)" << std::endl;
0161             gSystem->Exit(1);
0162           }
0163           if (verbosity > 0)
0164           {
0165             std::cout << std::endl
0166                       << std::endl
0167                       << "\tID number " << imagnet << std::endl;
0168             std::cout << "magname: " << magname << std::endl;
0169             std::cout << "x: " << x << std::endl;
0170             std::cout << "y: " << y << std::endl;
0171             std::cout << "z: " << z << std::endl;
0172             std::cout << "inner_radius_zin: " << inner_radius_zin << std::endl;
0173             std::cout << "inner_radius_zout: " << inner_radius_zout << std::endl;
0174             std::cout << "outer_magnet_diameter: " << outer_magnet_diameter << std::endl;
0175             std::cout << "length: " << length << std::endl;
0176             std::cout << "angle: " << angle << std::endl;
0177             std::cout << "dipole_field_x: " << dipole_field_x << std::endl;
0178             std::cout << "fieldgradient: " << fieldgradient << std::endl;
0179           }
0180           if (!magname.compare(0, 1, "B"))
0181           {
0182             magtype = "DIPOLE";
0183           }
0184           else if (!magname.compare(0, 1, "Q"))
0185           {
0186             magtype = "QUADRUPOLE";
0187           }
0188           else if (!magname.compare(0, 1, "S"))
0189           {
0190             magtype = "SEXTUPOLE";
0191           }
0192           else
0193           {
0194             std::cout << "cannot decode magnet name " << magname << std::endl;
0195             gSystem->Exit(1);
0196           }
0197           // convert to our units (cm, deg)
0198           x *= 100.;
0199           y *= 100.;
0200           z *= 100.;
0201           length *= 100.;
0202           inner_radius_zin *= 100.;
0203           outer_magnet_diameter *= 100.;
0204           angle = (angle / M_PI * 180.) / 1000.;  // given in mrad
0205
0206           //------------------------
0207
0208           if (magnetlist.empty() || magnetlist.contains(imagnet))
0209           {
0210             bl = new BeamLineMagnetSubsystem("BEAMLINEMAGNET", imagnet);
0211             bl->set_double_param("field_y", MagFieldFlip(dipole_field_x));
0212             bl->set_double_param("fieldgradient", MagFieldFlip(fieldgradient));
0213             bl->set_string_param("magtype", magtype);
0214             bl->set_double_param("length", length);
0215             bl->set_double_param("place_x", PosFlip(x));  // relative position to mother vol.
0216             bl->set_double_param("place_y", y);           // relative position to mother vol.
0217             if (z > 0)
0218             {
0219               bl->set_double_param("place_z", z - G4BEAMLINE::enclosure_center);  // relative position to mother vol.
0220             }
0221             else
0222             {
0223               bl->set_double_param("place_z", z + G4BEAMLINE::enclosure_center);  // relative position to mother vol.
0224             }
0225             bl->set_double_param("field_global_position_x", PosFlip(x));  // abs. position to world for field manager
0226             bl->set_double_param("field_global_position_y", y);           // abs. position to world for field manager
0227             bl->set_double_param("field_global_position_z", z);           // abs. position to world for field manager
0228             bl->set_double_param("rot_y", AngleFlip(angle));
0229             bl->set_double_param("field_global_rot_y", AngleFlip(angle));  // abs. rotation to world for field manager
0230             bl->set_double_param("inner_radius", inner_radius_zin);
0231             bl->set_double_param("outer_radius", outer_magnet_diameter / 2.);
0232             bl->set_double_param("skin_thickness", G4BEAMLINE::skin_thickness);
0233             bl->SetActive(AbsorberActive);
0234             bl->SetAbsorberActive(AbsorberActive);
0235             if (Enable::BEAMLINE_BLACKHOLE) bl->BlackHole();  // turn magnets into black holes
0236             if (z > 0)
0237             {
0238               bl->SetMotherSubsystem(G4BEAMLINE::ForwardBeamLineEnclosure);
0239             }
0240             else
0241             {
0242               bl->SetMotherSubsystem(G4BEAMLINE::BackwardBeamLineEnclosure);
0243             }
0244             bl->OverlapCheck(overlapCheck);
0245             bl->SuperDetector("BEAMLINEMAGNET");
0246             bl->Verbosity(verbosity);
0247             g4Reco->registerSubsystem(bl);
0248           }
0249           imagnet++;
0250           biggest_z = std::max(std::abs(z) + length, biggest_z);
0251         }
0252       }
0253     }
0254     infile.close();
0255   }
0256 }
0257
0258 void BeamLineDefineBeamPipe(PHG4Reco *g4Reco)
0259 {
0260   bool OverlapCheck = Enable::OVERLAPCHECK || Enable::BEAMLINE_OVERLAPCHECK;
0261   bool AbsorberActive = Enable::ABSORBER || Enable::BEAMLINE_ABSORBER;
0262   int verbosity = std::max(Enable::VERBOSITY, Enable::BEAMLINE_VERBOSITY);
0263
0264   const int ntube = 10;
0265   const std::string nm[ntube] = {"B00", "B01", "B10", "B11", "B20", "B21", "B30", "B31", "B32", "B33"};
0266   const double qlen[ntube] = {233.8, 233.8, 118.7, 118.7, 217.16, 217.16, 182.5, 182.5, 182.5, 182.5};
0267   const double qir[ntube] = {6.08, 6.08, 14.60, 14.60, 20.0, 20.0, 6.07, 6.07, 6.07, 6.07};
0268   const double qor[ntube] = {6.35, 6.35, 15.24, 15.24, 20.96, 20.96, 6.35, 6.35, 6.35, 6.35};
0269   const double qrot[ntube] = {0, 0, 0, 0, 0, 0, 1.074, -1.074, -1.074, 1.074};  // degree
0270   const double qxC[ntube] = {0, 0, 0, 0, 0, 0, 12.820, -12.820, 12.820, -12.820};
0271   const double qyC[ntube] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
0272   const double qzC[ntube] = {863.1, -863.1, 1474.470, -1474.470, 1642.4, -1642.4, 1843.2, 1843.2, -1843.2, -1843.2};
0273   for (int i = 0; i < ntube; i++)
0274   {
0275     std::string name = "beamPipe" + nm[i];
0276     PHG4CylinderSubsystem *pipe = new PHG4CylinderSubsystem(name, G4BEAMLINE::pipe_id_offset + i);
0277     if (Enable::BEAMLINE_BLACKHOLE) pipe->BlackHole();
0278     pipe->set_color(1, 0, 0, 1.);
0279     pipe->set_double_param("radius", qir[i]);
0280     pipe->set_double_param("thickness", qor[i] - qir[i]);
0281     pipe->set_double_param("length", qlen[i]);
0282     pipe->set_double_param("rot_y", qrot[i]);
0283     pipe->set_string_param("material", "G4_STAINLESS-STEEL");
0284     pipe->set_double_param("place_x", PosFlip(qxC[i]));
0285     pipe->set_double_param("place_y", qyC[i]);
0286     if (qzC[i] > 0)
0287     {
0288       pipe->set_double_param("place_z", qzC[i] - G4BEAMLINE::enclosure_center);  // relative position to mother vol.
0289     }
0290     else
0291     {
0292       pipe->set_double_param("place_z", qzC[i] + G4BEAMLINE::enclosure_center);  // relative position to mother vol.
0293     }
0294     if (AbsorberActive) pipe->SetActive();
0295     pipe->SuperDetector("PIPE");
0296     if (qzC[i] > 0)
0297     {
0298       pipe->SetMotherSubsystem(G4BEAMLINE::ForwardBeamLineEnclosure);
0299     }
0300     else
0301     {
0302       pipe->SetMotherSubsystem(G4BEAMLINE::BackwardBeamLineEnclosure);
0303     }
0304
0305     pipe->OverlapCheck(OverlapCheck);
0306     pipe->Verbosity(verbosity);
0307     g4Reco->registerSubsystem(pipe);
0308   }
0309
0310   // Roman Pot pipe
0311   const int nSec = 2;
0312   const double len[nSec] = {20.87, 20.87};
0313   const double ir1[nSec] = {7.14, 14.60};
0314   const double or1[nSec] = {7.77, 15.24};
0315   const double ir2[nSec] = {14.60, 7.14};
0316   const double or2[nSec] = {15.24, 7.77};
0317   const double xC[nSec] = {0, 0};
0318   const double yC[nSec] = {0, 0};
0319   const double zC[nSec] = {1394.25, -1394.25};
0320   for (int i = 0; i < nSec; i++)
0321   {
0322     std::string name = "beamPipeRP" + std::to_string(i);
0323     PHG4ConeSubsystem *pipe = new PHG4ConeSubsystem(name, G4BEAMLINE::roman_pot_pipe_id_offset + i);
0324     pipe->set_string_param("material", "G4_STAINLESS-STEEL");
0325     pipe->set_double_param("place_x", PosFlip(xC[i]));
0326     pipe->set_double_param("place_y", yC[i]);
0327     if (zC[i] > 0)
0328     {
0329       pipe->set_double_param("place_z", zC[i] - G4BEAMLINE::enclosure_center);
0330     }
0331     else
0332     {
0333       pipe->set_double_param("place_z", zC[i] + G4BEAMLINE::enclosure_center);
0334     }
0335     pipe->set_double_param("length", len[i]);
0336     pipe->set_double_param("rmin1", ir1[i]);
0337     pipe->set_double_param("rmin2", ir2[i]);
0338     pipe->set_double_param("rmax1", or1[i]);
0339     pipe->set_double_param("rmax2", or2[i]);
0340     if (zC[i] > 0)
0341     {
0342       pipe->SetMotherSubsystem(G4BEAMLINE::ForwardBeamLineEnclosure);
0343     }
0344     else
0345     {
0346       pipe->SetMotherSubsystem(G4BEAMLINE::BackwardBeamLineEnclosure);
0347     }
0348     if (AbsorberActive)
0349     {
0350       pipe->SetActive();
0351     }
0352     if (Enable::BEAMLINE_BLACKHOLE)
0353     {
0354       pipe->BlackHole();
0355     }
0356     pipe->set_color(1, 0, 0, 1.);
0357     pipe->SuperDetector("PIPE");
0358     pipe->OverlapCheck(OverlapCheck);
0359     pipe->Verbosity(verbosity);
0360     g4Reco->registerSubsystem(pipe);
0361   }
0362 }
0363
0364 float PosFlip(float pos)
0365 {
0366   return pos;
0367 };
0368 float AngleFlip(float angle)
0369 {
0370   return angle;
0371 };
0372 float MagFieldFlip(float Bfield)
0373 {
0374   return Bfield;
0375 };
0376
0377 #endif