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