g4simulation/g4detectors/PHG4OuterHcalField.cc

0001 // $Id: $
0002
0003 /*!
0004  * \file PHG4OuterHcalField.cc
0005  * \brief
0006  * \author Jin Huang <jhuang@bnl.gov>
0007  * \version $Revision:   $
0008  * \date $Date: $
0009  */
0010
0011 #include "PHG4OuterHcalField.h"
0012
0013 #include <TSystem.h>
0014
0015 #include <Geant4/G4Field.hh>  // for G4Field
0016 #include <Geant4/G4FieldManager.hh>
0017 #include <Geant4/G4PhysicalConstants.hh>
0018 #include <Geant4/G4SystemOfUnits.hh>
0019 #include <Geant4/G4TransportationManager.hh>
0020 #include <Geant4/G4Types.hh>  // for G4double, G4int
0021 #include <Geant4/G4Vector3D.hh>
0022
0023 #pragma GCC diagnostic push
0024 #pragma GCC diagnostic ignored "-Wshadow"
0025 #include <boost/stacktrace.hpp>
0026 #pragma GCC diagnostic pop
0027
0028 #include <cassert>  // for assert
0029 #include <cmath>    // for atan2, cos, sin, sqrt
0030 #include <cstdlib>  // for exit
0031 #include <iostream>
0032
0033 PHG4OuterHcalField::PHG4OuterHcalField(bool isInIron, G4int steelPlates,
0034                                        G4double scintiGap, G4double tiltAngle)
0035   : is_in_iron(isInIron)
0036   , n_steel_plates(steelPlates)
0037   , scinti_gap(scintiGap)
0038   , tilt_angle(tiltAngle)
0039 {
0040 }
0041
0042 void PHG4OuterHcalField::GetFieldValue(const double Point[4], double* Bfield) const
0043 {
0044   G4FieldManager* field_manager =
0045       G4TransportationManager::GetTransportationManager()->GetFieldManager();
0046
0047   if (!field_manager)
0048   {
0049     std::cout << "PHG4OuterHcalField::GetFieldValue"
0050               << " - Error! can not find field manager in G4TransportationManager"
0051               << std::endl;
0052     gSystem->Exit(1);
0053     exit(1);
0054   }
0055
0056   const G4Field* default_field = field_manager->GetDetectorField();
0057
0058   if (default_field)
0059   {
0060     default_field->GetFieldValue(Point, Bfield);
0061     // scale_factor for field component along the plate surface
0062     double x = Point[0];
0063     double y = Point[1];
0064     //      double z = Point[2];
0065
0066     assert(cos(tilt_angle) > 0);
0067     assert(n_steel_plates > 0);
0068
0069     // input 2D magnetic field vector
0070     const G4Vector3D B0(Bfield[0], Bfield[1], Bfield[2]);
0071     const G4Vector3D B0XY(Bfield[0], Bfield[1], 0);
0072     const G4Vector3D B0Z(0, 0, Bfield[2]);
0073
0074     const double R = sqrt(x * x + y * y);
0075     const double layer_RdPhi = R * twopi / n_steel_plates;
0076     const double layer_width = layer_RdPhi * cos(tilt_angle);
0077     const double gap_width = scinti_gap;
0078     if (gap_width >= layer_width)
0079     {
0080       std::cout << "PHG4OuterHcalField::GetFieldValue gap_width " << gap_width
0081                 << " < layer_width: " << layer_width
0082                 << " existing now, here is some debug info" << std::endl;
0083       std::cout << "coordinates: x: " << Point[0] / cm
0084                 << ", y: " << Point[1] / cm
0085                 << ", z: " << Point[2] / cm << std::endl;
0086       std::cout << "n_steel_plates: " << n_steel_plates << std::endl;
0087       std::cout << "Radius: " << R << std::endl;
0088       std::cout << "layer_RdPhi: " << layer_RdPhi << std::endl;
0089       std::cout << "layer_width: " << layer_width << std::endl;
0090       std::cout << "tilt_angle: " << tilt_angle << std::endl;
0091       std::cout << "And this is who called it:" << std::endl;
0092       std::cout << boost::stacktrace::stacktrace();
0093       std::cout << std::endl;
0094       return;
0095     }
0096     // sign definition of tilt_angle is rotation around the -z axis
0097     const G4Vector3D absorber_dir(cos(atan2(y, x) - tilt_angle),
0098                                   sin(atan2(y, x) - tilt_angle), 0);
0099     const G4Vector3D radial_dir(cos(atan2(y, x)), sin(atan2(y, x)), 0);
0100
0101     const double radial_flux_per_layer = layer_RdPhi * (B0XY.dot(radial_dir));
0102     double B_XY_mag = radial_flux_per_layer / (relative_permeability_absorber * (layer_width - gap_width) + relative_permeability_gap * gap_width);
0103     B_XY_mag *=
0104         is_in_iron ? relative_permeability_absorber : relative_permeability_gap;
0105
0106     const G4Vector3D B_New_XY = B_XY_mag * absorber_dir;
0107
0108     const double z_flux_per_layer = layer_width * B0Z.z();
0109     double B_Z_mag = z_flux_per_layer / (relative_permeability_absorber * (layer_width - gap_width) + relative_permeability_gap * gap_width);
0110     B_Z_mag *=
0111         is_in_iron ? relative_permeability_absorber : relative_permeability_gap;
0112     const G4Vector3D B_New_Z(0, 0, B_Z_mag);
0113
0114     // scale B_T component
0115     G4Vector3D B_New = B_New_Z + B_New_XY;
0116     Bfield[0] = B_New.x();
0117     Bfield[1] = B_New.y();
0118     Bfield[2] = B_New.z();
0119   }
0120   else
0121   {
0122     std::cout << "PHG4OuterHcalField::GetFieldValue"
0123               << " - Error! can not find detecor field in field manager!"
0124               << std::endl;
0125     gSystem->Exit(1);
0126     exit(1);
0127   }
0128 }