ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4BeamlineMagnetDetector.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4BeamlineMagnetDetector.cc
2 
3 #include <phparameter/PHParameters.h>
4 
5 #include <g4main/PHG4Detector.h> // for PHG4Detector
6 #include <g4main/PHG4Utils.h>
7 
8 #include <phool/phool.h>
9 
10 #include <Geant4/G4ChordFinder.hh>
11 #include <Geant4/G4ClassicalRK4.hh>
12 #include <Geant4/G4FieldManager.hh>
13 #include <Geant4/G4LogicalVolume.hh>
14 #include <Geant4/G4Mag_UsualEqRhs.hh>
15 #include <Geant4/G4MagneticField.hh>
16 #include <Geant4/G4PVPlacement.hh>
17 #include <Geant4/G4PhysicalConstants.hh>
18 #include <Geant4/G4QuadrupoleMagField.hh>
19 #include <Geant4/G4RotationMatrix.hh>
20 #include <Geant4/G4String.hh> // for G4String
21 #include <Geant4/G4SystemOfUnits.hh>
22 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
23 #include <Geant4/G4Transform3D.hh> // for G4Transform3D
24 #include <Geant4/G4Tubs.hh>
25 #include <Geant4/G4Types.hh> // for G4double, G4bool
26 #include <Geant4/G4UniformMagField.hh>
27 #include <Geant4/G4VisAttributes.hh>
28 
29 #include <CLHEP/Units/SystemOfUnits.h> // for cm, deg, tesla, twopi, meter
30 
31 #include <cstdlib> // for exit
32 #include <iostream> // for operator<<, basic_ostream
33 
34 class G4Material;
35 class G4VSolid;
36 class PHCompositeNode;
37 class PHG4Subsystem;
38 
39 using namespace std;
40 
41 //_______________________________________________________________
43  : PHG4Detector(subsys, Node, dnam)
44  , params(parameters)
45  , magnet_physi(nullptr)
46  , cylinder_physi(nullptr)
47  , layer(lyr)
48 {
49 }
50 
51 //_______________________________________________________________
53 {
54  if (volume == magnet_physi)
55  {
56  return true;
57  }
58  return false;
59 }
60 
61 //_______________________________________________________________
63 {
64  G4Material *TrackerMaterial = GetDetectorMaterial(params->get_string_param("material"));
65 
66  G4VisAttributes *fieldVis = new G4VisAttributes();
67  PHG4Utils::SetColour(fieldVis, "BlackHole");
68  fieldVis->SetVisibility(false);
69  fieldVis->SetForceSolid(false);
70 
71  G4VisAttributes *siliconVis = new G4VisAttributes();
72  if (params->get_int_param("blackhole"))
73  {
74  PHG4Utils::SetColour(siliconVis, "BlackHole");
75  siliconVis->SetVisibility(false);
76  siliconVis->SetForceSolid(false);
77  }
78  else
79  {
80  PHG4Utils::SetColour(siliconVis, params->get_string_param("material"));
81  siliconVis->SetVisibility(true);
82  siliconVis->SetForceSolid(true);
83  }
84 
85  /* Define origin vector (center of magnet) */
87  params->get_double_param("place_y") * cm,
88  params->get_double_param("place_z") * cm);
89 
90  /* Define magnet rotation matrix */
91  G4RotationMatrix *rotm = new G4RotationMatrix();
92  rotm->rotateX(params->get_double_param("rot_x") * deg);
93  rotm->rotateY(params->get_double_param("rot_y") * deg);
94  rotm->rotateZ(params->get_double_param("rot_z") * deg);
95 
96  /* Creating a magnetic field */
97  G4MagneticField *magField = nullptr;
98 
99  string magnettype = params->get_string_param("magtype");
100  if (magnettype == "dipole")
101  {
102  G4double fieldValue = params->get_double_param("field_y") * tesla;
103  magField = new G4UniformMagField(G4ThreeVector(0., fieldValue, 0.));
104 
105  if (Verbosity() > 0)
106  cout << "Creating DIPOLE with field " << fieldValue << " and name " << GetName() << endl;
107  }
108  else if (magnettype == "quadrupole")
109  {
110  G4double fieldGradient = params->get_double_param("fieldgradient") * tesla / meter;
111 
112  /* G4MagneticField::GetFieldValue( pos*, B* ) uses GLOBAL coordinates, not local.
113  * Therefore, place magnetic field center at the correct location and angle for the
114  * magnet AND do the same transformations for the logical volume (see below). */
115  magField = new G4QuadrupoleMagField(fieldGradient, origin, rotm);
116  // magField = new PHG4QuadrupoleMagField ( fieldGradient, origin, rotm );
117 
118  if (Verbosity() > 0)
119  {
120  cout << "Creating QUADRUPOLE with gradient " << fieldGradient << " and name " << GetName() << endl;
121  cout << "at x, y, z = " << origin.x() << " , " << origin.y() << " , " << origin.z() << endl;
122  cout << "with rotation around x, y, z axis by: " << rotm->phiX() << ", " << rotm->phiY() << ", " << rotm->phiZ() << endl;
123  }
124  }
125 
126  if (!magField)
127  {
128  cout << PHWHERE << " Aborting, no magnetic field specified for " << GetName() << endl;
129  exit(1);
130  }
131 
132  /* Set up Geant4 field manager */
133  G4Mag_UsualEqRhs *localEquation = new G4Mag_UsualEqRhs(magField);
134  G4ClassicalRK4 *localStepper = new G4ClassicalRK4(localEquation);
135  G4double minStep = 0.25 * mm; // minimal step, 1 mm is default
136  G4ChordFinder *localChordFinder = new G4ChordFinder(magField, minStep, localStepper);
137 
138  G4FieldManager *fieldMgr = new G4FieldManager();
139  fieldMgr->SetDetectorField(magField);
140  fieldMgr->SetChordFinder(localChordFinder);
141 
142  /* Add volume with magnetic field */
143  double radius = params->get_double_param("radius") * cm;
144  double thickness = params->get_double_param("thickness") * cm;
145 
146  G4VSolid *magnet_solid = new G4Tubs(GetName(),
147  0,
148  radius + thickness,
149  params->get_double_param("length") * cm / 2., 0, twopi);
150 
151  G4LogicalVolume *magnet_logic = new G4LogicalVolume(magnet_solid,
152  GetDetectorMaterial("G4_Galactic"),
153  GetName(),
154  0, 0, 0);
155  magnet_logic->SetVisAttributes(fieldVis);
156 
157  /* Set field manager for logical volume */
158  G4bool allLocal = true;
159  magnet_logic->SetFieldManager(fieldMgr, allLocal);
160 
161  /* create magnet physical volume */
163  G4ThreeVector(params->get_double_param("place_x") * cm,
164  params->get_double_param("place_y") * cm,
165  params->get_double_param("place_z") * cm)),
166  magnet_logic,
167  GetName(),
168  logicMother, 0, false, OverlapCheck());
169 
170  /* Add volume with solid magnet material */
171  G4VSolid *cylinder_solid = new G4Tubs(G4String(GetName().append("_Solid")),
172  radius,
173  radius + thickness,
174  params->get_double_param("length") * cm / 2., 0, twopi);
175  G4LogicalVolume *cylinder_logic = new G4LogicalVolume(cylinder_solid,
176  TrackerMaterial,
177  G4String(GetName()),
178  0, 0, 0);
179  cylinder_logic->SetVisAttributes(siliconVis);
180 
181  cylinder_physi = new G4PVPlacement(0, G4ThreeVector(0, 0, 0),
182  cylinder_logic,
183  G4String(GetName().append("_Solid")),
184  magnet_logic, 0, false, OverlapCheck());
185 }