g4simulation/g4detectors/BeamLineMagnetDetector.cc

0001 #include "BeamLineMagnetDetector.h"
0002 #include "BeamLineMagnetDisplayAction.h"
0003
0004 #include <phparameter/PHParameters.h>
0005
0006 #include <g4main/PHG4Detector.h>       // for PHG4Detector
0007 #include <g4main/PHG4DisplayAction.h>  // for PHG4DisplayAction
0008 #include <g4main/PHG4Subsystem.h>
0009
0010 #include <phool/phool.h>
0011
0012 #include <TSystem.h>
0013
0014 #include <Geant4/G4ChordFinder.hh>
0015 #include <Geant4/G4ClassicalRK4.hh>
0016 #include <Geant4/G4FieldManager.hh>
0017 #include <Geant4/G4LogicalVolume.hh>
0018 #include <Geant4/G4Mag_UsualEqRhs.hh>
0019 #include <Geant4/G4MagneticField.hh>
0020 #include <Geant4/G4PVPlacement.hh>
0021 #include <Geant4/G4PhysicalConstants.hh>
0022 #include <Geant4/G4QuadrupoleMagField.hh>
0023 #include <Geant4/G4RotationMatrix.hh>
0024 #include <Geant4/G4SystemOfUnits.hh>
0025 #include <Geant4/G4ThreeVector.hh>  // for G4ThreeVector
0026 #include <Geant4/G4Transform3D.hh>  // for G4Transform3D
0027 #include <Geant4/G4Tubs.hh>
0028 #include <Geant4/G4Types.hh>  // for G4double, G4bool
0029 #include <Geant4/G4UniformMagField.hh>
0030
0031 #include <CLHEP/Units/SystemOfUnits.h>  // for cm, deg, tesla, twopi, meter
0032
0033 #include <cassert>
0034 #include <cstdlib>   // for exit
0035 #include <iostream>  // for operator<<, basic_ostream
0036
0037 class G4VSolid;
0038 class PHCompositeNode;
0039
0040 //_______________________________________________________________
0041 BeamLineMagnetDetector::BeamLineMagnetDetector(PHG4Subsystem *subsys, PHCompositeNode *Node, PHParameters *parameters, const std::string &dnam, const int magnet_id)
0042   : PHG4Detector(subsys, Node, dnam)
0043   , m_Params(parameters)
0044   , m_DisplayAction(dynamic_cast<BeamLineMagnetDisplayAction *>(subsys->GetDisplayAction()))
0045   , m_MagnetId(magnet_id)
0046 {
0047 }
0048
0049 //_______________________________________________________________
0050 int BeamLineMagnetDetector::IsInBeamLineMagnet(const G4VPhysicalVolume *volume) const
0051 {
0052   if (volume == magnet_physi)
0053   {
0054     return 1;
0055   }
0056   if (volume == magnet_core_physi)
0057   {
0058     return -2;
0059   }
0060   if (volume == magnet_iron_physi)
0061   {
0062     return -1;
0063   }
0064   return 0;
0065 }
0066
0067 //_______________________________________________________________
0068 void BeamLineMagnetDetector::ConstructMe(G4LogicalVolume *logicMother)
0069 {
0070   /* Define origin vector (center of magnet) */
0071   G4ThreeVector origin(m_Params->get_double_param("place_x") * cm,
0072                        m_Params->get_double_param("place_y") * cm,
0073                        m_Params->get_double_param("place_z") * cm);
0074
0075   /* Define magnet rotation matrix */
0076   G4RotationMatrix *rotm = new G4RotationMatrix();
0077   rotm->rotateX(m_Params->get_double_param("rot_x") * deg);
0078   rotm->rotateY(m_Params->get_double_param("rot_y") * deg);
0079   rotm->rotateZ(m_Params->get_double_param("rot_z") * deg);
0080
0081   // creating mother volume (cylinder for the time being)
0082   G4VSolid *magnet_mother_solid = new G4Tubs(GetName().append("_Mother"),
0083                                              0,
0084                                              m_Params->get_double_param("outer_radius") * cm,
0085                                              m_Params->get_double_param("length") * cm / 2., 0, twopi);
0086   G4LogicalVolume *magnet_mother_logic = new G4LogicalVolume(magnet_mother_solid,
0087                                                              GetDetectorMaterial("G4_Galactic"),
0088                                                              GetName(),
0089                                                              nullptr, nullptr, nullptr);
0090
0091   m_DisplayAction->AddVolume(magnet_mother_logic, "OFF");
0092   new G4PVPlacement(G4Transform3D(*rotm, origin),
0093                     magnet_mother_logic,
0094                     GetName().append("_Mother"),
0095                     logicMother, false, m_MagnetId, OverlapCheck());
0096
0097   G4ThreeVector field_origin(origin);
0098   G4RotationMatrix *field_rotm = rotm;
0099
0100   // mother subsys ! = the world and therefore need to explicitly assign global geometry for field mamnagers
0101   if (GetMySubsystem()->GetMotherSubsystem())
0102   {
0103     if (Verbosity() > 0)
0104     {
0105       std::cout << __PRETTY_FUNCTION__ << ": set field using the global coordinate system, as the magnet is a daughter vol. of "
0106                 << GetMySubsystem()->GetMotherSubsystem()->Name() << std::endl;
0107     }
0108     /* Define origin vector (center of magnet) */
0109     // abs. position to world for field manager
0110     field_origin = G4ThreeVector(m_Params->get_double_param("field_global_position_x") * cm,
0111                                  m_Params->get_double_param("field_global_position_y") * cm,
0112                                  m_Params->get_double_param("field_global_position_z") * cm);
0113
0114     /* Define magnet rotation matrix */
0115     // abs. rotation to world for field manager
0116     field_rotm = new G4RotationMatrix();
0117     rotm->rotateX(m_Params->get_double_param("field_global_rot_x") * deg);
0118     rotm->rotateY(m_Params->get_double_param("field_global_rot_y") * deg);
0119     rotm->rotateZ(m_Params->get_double_param("field_global_rot_z") * deg);
0120   }
0121
0122   std::string magnettype = m_Params->get_string_param("magtype");
0123   if (magnettype == "DIPOLE")
0124   {
0125     G4ThreeVector field(m_Params->get_double_param("field_x") * tesla,
0126                         m_Params->get_double_param("field_y") * tesla,
0127                         m_Params->get_double_param("field_z") * tesla);
0128     // magnets can be rotated in y
0129     //     field.rotateY(m_Params->get_double_param("rot_y") * deg);
0130     // apply consistent geometry ops to field and detectors
0131     field.transform(*field_rotm);
0132     m_magField = new G4UniformMagField(field);
0133     if (Verbosity() > 0)
0134     {
0135       std::cout << "Creating DIPOLE with field x: " << field.x() / tesla
0136                 << ", y: " << field.y() / tesla
0137                 << ", z: " << field.z() / tesla
0138                 << " and name " << GetName() << std::endl;
0139     }
0140   }
0141   else if (magnettype == "QUADRUPOLE")
0142   {
0143     G4double fieldGradient = m_Params->get_double_param("fieldgradient") * tesla / meter;
0144
0145     /* G4MagneticField::GetFieldValue( pos*, B* ) uses GLOBAL coordinates, not local.
0146      * Therefore, place magnetic field center at the correct location and angle for the
0147      * magnet AND do the same transformations for the logical volume (see below). */
0148     m_magField = new G4QuadrupoleMagField(fieldGradient, field_origin, field_rotm);
0149
0150     if (Verbosity() > 0)
0151     {
0152       std::cout << "Creating QUADRUPOLE with gradient " << fieldGradient << " and name " << GetName() << std::endl;
0153       std::cout << "at x, y, z = " << field_origin.x() << " , " << field_origin.y() << " , " << field_origin.z() << std::endl;
0154       std::cout << "with rotation around x, y, z axis  by: " << field_rotm->phiX() << ", " << field_rotm->phiY() << ", " << field_rotm->phiZ() << std::endl;
0155     }
0156   }
0157   else
0158   {
0159     std::cout << __PRETTY_FUNCTION__ << " : Fatal error, magnet type of " << magnettype << " is not yet supported" << std::endl;
0160     exit(1);
0161   }
0162
0163   if (!m_magField && Verbosity() > 0)
0164   {
0165     std::cout << PHWHERE << " No magnetic field specified for " << GetName()
0166               << " of type " << magnettype << std::endl;
0167   }
0168
0169   /* Add volume with solid magnet material */
0170   G4VSolid *magnet_iron_solid = new G4Tubs(GetName().append("_Solid"),
0171                                            m_Params->get_double_param("inner_radius") * cm,
0172                                            m_Params->get_double_param("outer_radius") * cm,
0173                                            m_Params->get_double_param("length") * cm / 2., 0, twopi);
0174   G4LogicalVolume *magnet_iron_logic = new G4LogicalVolume(magnet_iron_solid,
0175                                                            GetDetectorMaterial("G4_Fe"),
0176                                                            GetName(),
0177                                                            nullptr, nullptr, nullptr);
0178   m_DisplayAction->AddVolume(magnet_iron_logic, magnettype);
0179
0180   if (m_Params->get_double_param("skin_thickness") > 0)
0181   {
0182     if (m_Params->get_double_param("inner_radius") * cm + m_Params->get_double_param("skin_thickness") * cm >
0183         m_Params->get_double_param("outer_radius") * cm - m_Params->get_double_param("skin_thickness") * cm)
0184     {
0185       std::cout << "Magnet skin thickness " << m_Params->get_double_param("skin_thickness") << " too large, exceeds 2xmagnet thickness "
0186                 << m_Params->get_double_param("outer_radius") * cm - m_Params->get_double_param("inner_radius") * cm
0187                 << std::endl;
0188       gSystem->Exit(1);
0189     }
0190
0191     G4VSolid *magnet_core_solid = new G4Tubs(GetName().append("_Core"),
0192                                              m_Params->get_double_param("inner_radius") * cm + m_Params->get_double_param("skin_thickness") * cm,
0193                                              m_Params->get_double_param("outer_radius") * cm - m_Params->get_double_param("skin_thickness") * cm,
0194                                              (m_Params->get_double_param("length") * cm - m_Params->get_double_param("skin_thickness") * cm) / 2., 0, twopi);
0195     G4LogicalVolume *magnet_core_logic = new G4LogicalVolume(magnet_core_solid,
0196                                                              GetDetectorMaterial("G4_Fe"),
0197                                                              GetName(),
0198                                                              nullptr, nullptr, nullptr);
0199     m_DisplayAction->AddVolume(magnet_core_logic, magnettype);
0200
0201     magnet_core_physi = new G4PVPlacement(nullptr, G4ThreeVector(0, 0, 0),
0202                                           magnet_core_logic,
0203                                           GetName().append("_Solid"),
0204                                           magnet_iron_logic, false, m_MagnetId, OverlapCheck());
0205   }
0206   magnet_iron_physi = new G4PVPlacement(nullptr, G4ThreeVector(0, 0, 0),
0207                                         magnet_iron_logic,
0208                                         GetName().append("_Solid"),
0209                                         magnet_mother_logic, false, m_MagnetId, OverlapCheck());
0210
0211   G4VSolid *magnet_field_solid = new G4Tubs(GetName().append("_Field_Solid"),
0212                                             0,
0213                                             m_Params->get_double_param("inner_radius") * cm,
0214                                             m_Params->get_double_param("length") * cm / 2., 0, twopi);
0215
0216   m_magnetFieldLogic = new G4LogicalVolume(magnet_field_solid,
0217                                            GetDetectorMaterial("G4_Galactic"),
0218                                            GetName(),
0219                                            nullptr, nullptr, nullptr);
0220
0221   // allow installing new detector inside the magnet, magnet_field_logic
0222   PHG4Subsystem *mysys = GetMySubsystem();
0223   mysys->SetLogicalVolume(m_magnetFieldLogic);
0224
0225   m_DisplayAction->AddVolume(m_magnetFieldLogic, "FIELDVOLUME");
0226
0227   /* create magnet physical volume */
0228
0229   magnet_physi = new G4PVPlacement(nullptr, G4ThreeVector(0, 0, 0),
0230                                    m_magnetFieldLogic,
0231                                    GetName(),
0232                                    magnet_mother_logic, false, m_MagnetId, OverlapCheck());
0233 }
0234
0235 //! Optional PostConstruction call after all geometry is constructed
0236 void BeamLineMagnetDetector::PostConstruction()
0237 {
0238   /* Set field manager for logical volume */
0239   // This has to be in PostConstruction() in order to apply to all daughter vol. including daughter subsystems
0240
0241   assert(m_magField);
0242   assert(m_magnetFieldLogic);
0243
0244   if (Verbosity() > 0)
0245   {
0246     std::cout << __PRETTY_FUNCTION__ << ": set field to vol " << m_magnetFieldLogic->GetName() << " that include "
0247               << m_magnetFieldLogic->GetNoDaughters() << " daughter vols." << std::endl;
0248   }
0249
0250   /* Set up Geant4 field manager */
0251   G4Mag_UsualEqRhs *localEquation = new G4Mag_UsualEqRhs(m_magField);
0252   G4ClassicalRK4 *localStepper = new G4ClassicalRK4(localEquation);
0253   G4double minStep = 0.25 * mm;  // minimal step, 1 mm is default
0254   G4ChordFinder *localChordFinder = new G4ChordFinder(m_magField, minStep, localStepper);
0255
0256   G4FieldManager *fieldMgr = new G4FieldManager();
0257   fieldMgr->SetDetectorField(m_magField);
0258   fieldMgr->SetChordFinder(localChordFinder);
0259   m_magnetFieldLogic->SetFieldManager(fieldMgr, true);
0260 }