Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BeamLineMagnetDetector.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file BeamLineMagnetDetector.cc
3 
4 #include <phparameter/PHParameters.h>
5 
6 #include <g4main/PHG4Detector.h> // for PHG4Detector
7 #include <g4main/PHG4DisplayAction.h> // for PHG4DisplayAction
8 #include <g4main/PHG4Subsystem.h>
9 
10 #include <phool/phool.h>
11 
12 #include <TSystem.h>
13 
14 #include <Geant4/G4ChordFinder.hh>
15 #include <Geant4/G4ClassicalRK4.hh>
16 #include <Geant4/G4FieldManager.hh>
17 #include <Geant4/G4LogicalVolume.hh>
18 #include <Geant4/G4Mag_UsualEqRhs.hh>
19 #include <Geant4/G4MagneticField.hh>
20 #include <Geant4/G4PVPlacement.hh>
21 #include <Geant4/G4PhysicalConstants.hh>
22 #include <Geant4/G4QuadrupoleMagField.hh>
23 #include <Geant4/G4RotationMatrix.hh>
24 #include <Geant4/G4SystemOfUnits.hh>
25 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
26 #include <Geant4/G4Transform3D.hh> // for G4Transform3D
27 #include <Geant4/G4Tubs.hh>
28 #include <Geant4/G4Types.hh> // for G4double, G4bool
29 #include <Geant4/G4UniformMagField.hh>
30 
31 #include <CLHEP/Units/SystemOfUnits.h> // for cm, deg, tesla, twopi, meter
32 
33 #include <cassert>
34 #include <cstdlib> // for exit
35 #include <iostream> // for operator<<, basic_ostream
36 
37 class G4VSolid;
38 class PHCompositeNode;
39 
40 //_______________________________________________________________
42  : PHG4Detector(subsys, Node, dnam)
43  , m_Params(parameters)
44  , m_DisplayAction(dynamic_cast<BeamLineMagnetDisplayAction *>(subsys->GetDisplayAction()))
45  , m_MagnetId(magnet_id)
46 {
47 }
48 
49 //_______________________________________________________________
51 {
52  if (volume == magnet_physi)
53  {
54  return 1;
55  }
56  if (volume == magnet_core_physi)
57  {
58  return -2;
59  }
60  if (volume == magnet_iron_physi)
61  {
62  return -1;
63  }
64  return 0;
65 }
66 
67 //_______________________________________________________________
68 void BeamLineMagnetDetector::ConstructMe(G4LogicalVolume *logicMother)
69 {
70  /* Define origin vector (center of magnet) */
71  G4ThreeVector origin(m_Params->get_double_param("place_x") * cm,
72  m_Params->get_double_param("place_y") * cm,
73  m_Params->get_double_param("place_z") * cm);
74 
75  /* Define magnet rotation matrix */
76  G4RotationMatrix *rotm = new G4RotationMatrix();
77  rotm->rotateX(m_Params->get_double_param("rot_x") * deg);
78  rotm->rotateY(m_Params->get_double_param("rot_y") * deg);
79  rotm->rotateZ(m_Params->get_double_param("rot_z") * deg);
80 
81  // creating mother volume (cylinder for the time being)
82  G4VSolid *magnet_mother_solid = new G4Tubs(GetName().append("_Mother"),
83  0,
84  m_Params->get_double_param("outer_radius") * cm,
85  m_Params->get_double_param("length") * cm / 2., 0, twopi);
86  G4LogicalVolume *magnet_mother_logic = new G4LogicalVolume(magnet_mother_solid,
87  GetDetectorMaterial("G4_Galactic"),
88  GetName(),
89  nullptr, nullptr, nullptr);
90 
91  m_DisplayAction->AddVolume(magnet_mother_logic, "OFF");
92  new G4PVPlacement(G4Transform3D(*rotm, origin),
93  magnet_mother_logic,
94  GetName().append("_Mother"),
95  logicMother, false, m_MagnetId, OverlapCheck());
96 
97  G4ThreeVector field_origin(origin);
98  G4RotationMatrix *field_rotm = rotm;
99 
100  // mother subsys ! = the world and therefore need to explicitly assign global geometry for field mamnagers
101  if (GetMySubsystem()->GetMotherSubsystem())
102  {
103  if (Verbosity() > 0)
104  {
105  std::cout << __PRETTY_FUNCTION__ << ": set field using the global coordinate system, as the magnet is a daughter vol. of "
106  << GetMySubsystem()->GetMotherSubsystem()->Name() << std::endl;
107  }
108  /* Define origin vector (center of magnet) */
109  // abs. position to world for field manager
110  field_origin = G4ThreeVector(m_Params->get_double_param("field_global_position_x") * cm,
111  m_Params->get_double_param("field_global_position_y") * cm,
112  m_Params->get_double_param("field_global_position_z") * cm);
113 
114  /* Define magnet rotation matrix */
115  // abs. rotation to world for field manager
116  field_rotm = new G4RotationMatrix();
117  rotm->rotateX(m_Params->get_double_param("field_global_rot_x") * deg);
118  rotm->rotateY(m_Params->get_double_param("field_global_rot_y") * deg);
119  rotm->rotateZ(m_Params->get_double_param("field_global_rot_z") * deg);
120  }
121 
122  std::string magnettype = m_Params->get_string_param("magtype");
123  if (magnettype == "DIPOLE")
124  {
125  G4ThreeVector field(m_Params->get_double_param("field_x") * tesla,
126  m_Params->get_double_param("field_y") * tesla,
127  m_Params->get_double_param("field_z") * tesla);
128  // magnets can be rotated in y
129  // field.rotateY(m_Params->get_double_param("rot_y") * deg);
130  // apply consistent geometry ops to field and detectors
131  field.transform(*field_rotm);
132  m_magField = new G4UniformMagField(field);
133  if (Verbosity() > 0)
134  {
135  std::cout << "Creating DIPOLE with field x: " << field.x() / tesla
136  << ", y: " << field.y() / tesla
137  << ", z: " << field.z() / tesla
138  << " and name " << GetName() << std::endl;
139  }
140  }
141  else if (magnettype == "QUADRUPOLE")
142  {
143  G4double fieldGradient = m_Params->get_double_param("fieldgradient") * tesla / meter;
144 
145  /* G4MagneticField::GetFieldValue( pos*, B* ) uses GLOBAL coordinates, not local.
146  * Therefore, place magnetic field center at the correct location and angle for the
147  * magnet AND do the same transformations for the logical volume (see below). */
148  m_magField = new G4QuadrupoleMagField(fieldGradient, field_origin, field_rotm);
149 
150  if (Verbosity() > 0)
151  {
152  std::cout << "Creating QUADRUPOLE with gradient " << fieldGradient << " and name " << GetName() << std::endl;
153  std::cout << "at x, y, z = " << field_origin.x() << " , " << field_origin.y() << " , " << field_origin.z() << std::endl;
154  std::cout << "with rotation around x, y, z axis by: " << field_rotm->phiX() << ", " << field_rotm->phiY() << ", " << field_rotm->phiZ() << std::endl;
155  }
156  }
157  else
158  {
159  std::cout << __PRETTY_FUNCTION__ << " : Fatal error, magnet type of " << magnettype << " is not yet supported" << std::endl;
160  exit(1);
161  }
162 
163  if (!m_magField && Verbosity() > 0)
164  {
165  std::cout << PHWHERE << " No magnetic field specified for " << GetName()
166  << " of type " << magnettype << std::endl;
167  }
168 
169  /* Add volume with solid magnet material */
170  G4VSolid *magnet_iron_solid = new G4Tubs(GetName().append("_Solid"),
171  m_Params->get_double_param("inner_radius") * cm,
172  m_Params->get_double_param("outer_radius") * cm,
173  m_Params->get_double_param("length") * cm / 2., 0, twopi);
174  G4LogicalVolume *magnet_iron_logic = new G4LogicalVolume(magnet_iron_solid,
175  GetDetectorMaterial("G4_Fe"),
176  GetName(),
177  nullptr, nullptr, nullptr);
178  m_DisplayAction->AddVolume(magnet_iron_logic, magnettype);
179 
180  if (m_Params->get_double_param("skin_thickness") > 0)
181  {
182  if (m_Params->get_double_param("inner_radius") * cm + m_Params->get_double_param("skin_thickness") * cm >
183  m_Params->get_double_param("outer_radius") * cm - m_Params->get_double_param("skin_thickness") * cm)
184  {
185  std::cout << "Magnet skin thickness " << m_Params->get_double_param("skin_thickness") << " too large, exceeds 2xmagnet thickness "
186  << m_Params->get_double_param("outer_radius") * cm - m_Params->get_double_param("inner_radius") * cm
187  << std::endl;
188  gSystem->Exit(1);
189  }
190 
191  G4VSolid *magnet_core_solid = new G4Tubs(GetName().append("_Core"),
192  m_Params->get_double_param("inner_radius") * cm + m_Params->get_double_param("skin_thickness") * cm,
193  m_Params->get_double_param("outer_radius") * cm - m_Params->get_double_param("skin_thickness") * cm,
194  (m_Params->get_double_param("length") * cm - m_Params->get_double_param("skin_thickness") * cm) / 2., 0, twopi);
195  G4LogicalVolume *magnet_core_logic = new G4LogicalVolume(magnet_core_solid,
196  GetDetectorMaterial("G4_Fe"),
197  GetName(),
198  nullptr, nullptr, nullptr);
199  m_DisplayAction->AddVolume(magnet_core_logic, magnettype);
200 
201  magnet_core_physi = new G4PVPlacement(nullptr, G4ThreeVector(0, 0, 0),
202  magnet_core_logic,
203  GetName().append("_Solid"),
204  magnet_iron_logic, false, m_MagnetId, OverlapCheck());
205  }
206  magnet_iron_physi = new G4PVPlacement(nullptr, G4ThreeVector(0, 0, 0),
207  magnet_iron_logic,
208  GetName().append("_Solid"),
209  magnet_mother_logic, false, m_MagnetId, OverlapCheck());
210 
211  G4VSolid *magnet_field_solid = new G4Tubs(GetName().append("_Field_Solid"),
212  0,
213  m_Params->get_double_param("inner_radius") * cm,
214  m_Params->get_double_param("length") * cm / 2., 0, twopi);
215 
216  m_magnetFieldLogic = new G4LogicalVolume(magnet_field_solid,
217  GetDetectorMaterial("G4_Galactic"),
218  GetName(),
219  nullptr, nullptr, nullptr);
220 
221  // allow installing new detector inside the magnet, magnet_field_logic
222  PHG4Subsystem *mysys = GetMySubsystem();
223  mysys->SetLogicalVolume(m_magnetFieldLogic);
224 
225  m_DisplayAction->AddVolume(m_magnetFieldLogic, "FIELDVOLUME");
226 
227  /* create magnet physical volume */
228 
229  magnet_physi = new G4PVPlacement(nullptr, G4ThreeVector(0, 0, 0),
230  m_magnetFieldLogic,
231  GetName(),
232  magnet_mother_logic, false, m_MagnetId, OverlapCheck());
233 }
234 
237 {
238  /* Set field manager for logical volume */
239  // This has to be in PostConstruction() in order to apply to all daughter vol. including daughter subsystems
240 
243 
244  if (Verbosity() > 0)
245  {
246  std::cout << __PRETTY_FUNCTION__ << ": set field to vol " << m_magnetFieldLogic->GetName() << " that include "
247  << m_magnetFieldLogic->GetNoDaughters() << " daughter vols." << std::endl;
248  }
249 
250  /* Set up Geant4 field manager */
251  G4Mag_UsualEqRhs *localEquation = new G4Mag_UsualEqRhs(m_magField);
252  G4ClassicalRK4 *localStepper = new G4ClassicalRK4(localEquation);
253  G4double minStep = 0.25 * mm; // minimal step, 1 mm is default
254  G4ChordFinder *localChordFinder = new G4ChordFinder(m_magField, minStep, localStepper);
255 
256  G4FieldManager *fieldMgr = new G4FieldManager();
257  fieldMgr->SetDetectorField(m_magField);
258  fieldMgr->SetChordFinder(localChordFinder);
259  m_magnetFieldLogic->SetFieldManager(fieldMgr, true);
260 }