ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4OuterHcalField.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4OuterHcalField.cc
1 // $Id: $
2 
11 #include "PHG4OuterHcalField.h"
12 
13 #include <TSystem.h>
14 
15 #include <Geant4/G4Field.hh> // for G4Field
16 #include <Geant4/G4FieldManager.hh>
17 #include <Geant4/G4PhysicalConstants.hh>
18 #include <Geant4/G4SystemOfUnits.hh>
19 #include <Geant4/G4TransportationManager.hh>
20 #include <Geant4/G4Types.hh> // for G4double, G4int
21 #include <Geant4/G4Vector3D.hh>
22 
23 #include <cassert> // for assert
24 #include <cmath> // for atan2, cos, sin, sqrt
25 #include <cstdlib> // for exit
26 #include <iostream>
27 
28 PHG4OuterHcalField::PHG4OuterHcalField(bool isInIron, G4int steelPlates,
29  G4double scintiGap, G4double tiltAngle)
30  : is_in_iron(isInIron)
31  , n_steel_plates(steelPlates)
32  , scinti_gap(scintiGap)
33  , tilt_angle(tiltAngle)
34 {
35 }
36 
37 void PHG4OuterHcalField::GetFieldValue(const double Point[4], double* Bfield) const
38 {
39  G4FieldManager* field_manager =
41 
42  if (!field_manager)
43  {
44  std::cout << "PHG4OuterHcalField::GetFieldValue"
45  << " - Error! can not find field manager in G4TransportationManager"
46  << std::endl;
47  gSystem->Exit(1);
48  exit(1);
49  }
50 
51  const G4Field* default_field = field_manager->GetDetectorField();
52 
53  if (default_field)
54  {
55  default_field->GetFieldValue(Point, Bfield);
56 
57  // scale_factor for field component along the plate surface
58  double x = Point[0];
59  double y = Point[1];
60  // double z = Point[2];
61 
62  assert(cos(tilt_angle) > 0);
63  assert(n_steel_plates > 0);
64 
65  // input 2D magnetic field vector
66  const G4Vector3D B0(Bfield[0], Bfield[1], Bfield[2]);
67  const G4Vector3D B0XY(Bfield[0], Bfield[1], 0);
68  const G4Vector3D B0Z(0, 0, Bfield[2]);
69 
70  const double R = sqrt(x * x + y * y);
71  const double layer_RdPhi = R * twopi / n_steel_plates;
72  const double layer_width = layer_RdPhi * cos(tilt_angle);
73  const double gap_width = scinti_gap;
74  if (gap_width >= layer_width)
75  {
76  std::cout << "PHG4OuterHcalField::GetFieldValue gap_width " << gap_width
77  << " < layer_width: " << layer_width
78  << " existing now, here is some debug info" << std::endl;
79  std::cout << "coordinates: x: " << Point[0] / cm
80  << ", y: " << Point[1] / cm
81  << ", z: " << Point[2] / cm << std::endl;
82  std::cout << "n_steel_plates: " << n_steel_plates << std::endl;
83  std::cout << "Radius: " << R << std::endl;
84  std::cout << "layer_RdPhi: " << layer_RdPhi << std::endl;
85  std::cout << "layer_width: " << layer_width << std::endl;
86  std::cout << "tilt_angle: " << tilt_angle << std::endl;
87  gSystem->Exit(1);
88  exit(1);
89  }
90  // sign definition of tilt_angle is rotation around the -z axis
91  const G4Vector3D absorber_dir(cos(atan2(y, x) - tilt_angle),
92  sin(atan2(y, x) - tilt_angle), 0);
93  const G4Vector3D radial_dir(cos(atan2(y, x)), sin(atan2(y, x)), 0);
94 
95  const double radial_flux_per_layer = layer_RdPhi * (B0XY.dot(radial_dir));
96  double B_XY_mag = radial_flux_per_layer / (relative_permeability_absorber * (layer_width - gap_width) + relative_permeability_gap * gap_width);
97  B_XY_mag *=
99 
100  const G4Vector3D B_New_XY = B_XY_mag * absorber_dir;
101 
102  const double z_flux_per_layer = layer_width * B0Z.z();
103  double B_Z_mag = z_flux_per_layer / (relative_permeability_absorber * (layer_width - gap_width) + relative_permeability_gap * gap_width);
104  B_Z_mag *=
106  const G4Vector3D B_New_Z(0, 0, B_Z_mag);
107 
108  // scale B_T component
109  G4Vector3D B_New = B_New_Z + B_New_XY;
110  Bfield[0] = B_New.x();
111  Bfield[1] = B_New.y();
112  Bfield[2] = B_New.z();
113  }
114  else
115  {
116  std::cout << "PHG4OuterHcalField::GetFieldValue"
117  << " - Error! can not find detecor field in field manager!"
118  << std::endl;
119  gSystem->Exit(1);
120  exit(1);
121  }
122 }