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