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 #include <boost/stacktrace.hpp>
0024
0025 #include <cassert>  // for assert
0026 #include <cmath>    // for atan2, cos, sin, sqrt
0027 #include <cstdlib>  // for exit
0028 #include <iostream>
0029
0030 PHG4OuterHcalField::PHG4OuterHcalField(bool isInIron, G4int steelPlates,
0031                                        G4double scintiGap, G4double tiltAngle)
0032   : is_in_iron(isInIron)
0033   , n_steel_plates(steelPlates)
0034   , scinti_gap(scintiGap)
0035   , tilt_angle(tiltAngle)
0036 {
0037 }
0038
0039 void PHG4OuterHcalField::GetFieldValue(const double Point[4], double* Bfield) const
0040 {
0041   G4FieldManager* field_manager = G4TransportationManager::GetTransportationManager()->GetFieldManager();
0042
0043   if (!field_manager)
0044   {
0045     std::cout << "PHG4OuterHcalField::GetFieldValue"
0046               << " - Error! can not find field manager in G4TransportationManager"
0047               << std::endl;
0048     gSystem->Exit(1);
0049     exit(1);
0050   }
0051
0052   const G4Field* default_field = field_manager->GetDetectorField();
0053
0054   if (default_field)
0055   {
0056     default_field->GetFieldValue(Point, Bfield);
0057     // scale_factor for field component along the plate surface
0058     double x = Point[0];
0059     double y = Point[1];
0060     //      double z = Point[2];
0061
0062     assert(cos(tilt_angle) > 0);
0063     assert(n_steel_plates > 0);
0064
0065     // input 2D magnetic field vector
0066     const G4Vector3D B0(Bfield[0], Bfield[1], Bfield[2]);
0067     const G4Vector3D B0XY(Bfield[0], Bfield[1], 0);
0068     const G4Vector3D B0Z(0, 0, Bfield[2]);
0069
0070     const double R = sqrt(x * x + y * y);
0071     const double layer_RdPhi = R * twopi / n_steel_plates;
0072     const double layer_width = layer_RdPhi * cos(tilt_angle);
0073     const double gap_width = scinti_gap;
0074     if (gap_width >= layer_width)
0075     {
0076       std::cout << "PHG4OuterHcalField::GetFieldValue gap_width " << gap_width
0077                 << " < layer_width: " << layer_width
0078                 << " existing now, here is some debug info" << std::endl;
0079       std::cout << "coordinates: x: " << Point[0] / cm
0080                 << ", y: " << Point[1] / cm
0081                 << ", z: " << Point[2] / cm << std::endl;
0082       std::cout << "n_steel_plates: " << n_steel_plates << std::endl;
0083       std::cout << "Radius: " << R << std::endl;
0084       std::cout << "layer_RdPhi: " << layer_RdPhi << std::endl;
0085       std::cout << "layer_width: " << layer_width << std::endl;
0086       std::cout << "tilt_angle: " << tilt_angle << std::endl;
0087       std::cout << "And this is who called it:" << std::endl;
0088       std::cout << boost::stacktrace::stacktrace();
0089       std::cout << std::endl;
0090       return;
0091     }
0092     // sign definition of tilt_angle is rotation around the -z axis
0093     const G4Vector3D absorber_dir(cos(atan2(y, x) - tilt_angle),
0094                                   sin(atan2(y, x) - tilt_angle), 0);
0095     const G4Vector3D radial_dir(cos(atan2(y, x)), sin(atan2(y, x)), 0);
0096
0097     const double radial_flux_per_layer = layer_RdPhi * (B0XY.dot(radial_dir));
0098     double B_XY_mag = radial_flux_per_layer / (relative_permeability_absorber * (layer_width - gap_width) + relative_permeability_gap * gap_width);
0099     B_XY_mag *=
0100         is_in_iron ? relative_permeability_absorber : relative_permeability_gap;
0101
0102     const G4Vector3D B_New_XY = B_XY_mag * absorber_dir;
0103
0104     const double z_flux_per_layer = layer_width * B0Z.z();
0105     double B_Z_mag = z_flux_per_layer / (relative_permeability_absorber * (layer_width - gap_width) + relative_permeability_gap * gap_width);
0106     B_Z_mag *=
0107         is_in_iron ? relative_permeability_absorber : relative_permeability_gap;
0108     const G4Vector3D B_New_Z(0, 0, B_Z_mag);
0109
0110     // scale B_T component
0111     G4Vector3D B_New = B_New_Z + B_New_XY;
0112     Bfield[0] = B_New.x();
0113     Bfield[1] = B_New.y();
0114     Bfield[2] = B_New.z();
0115   }
0116   else
0117   {
0118     std::cout << "PHG4OuterHcalField::GetFieldValue"
0119               << " - Error! can not find detecor field in field manager!"
0120               << std::endl;
0121     gSystem->Exit(1);
0122     exit(1);
0123   }
0124 }