g4simulation/g4detectors/PHG4BeamlineMagnetDetector.cc

0001 #include "PHG4BeamlineMagnetDetector.h"
0002
0003 #include <phparameter/PHParameters.h>
0004
0005 #include <g4main/PHG4Detector.h>  // for PHG4Detector
0006 #include <g4main/PHG4Utils.h>
0007
0008 #include <phool/phool.h>
0009
0010 #include <Geant4/G4ChordFinder.hh>
0011 #include <Geant4/G4ClassicalRK4.hh>
0012 #include <Geant4/G4FieldManager.hh>
0013 #include <Geant4/G4LogicalVolume.hh>
0014 #include <Geant4/G4Mag_UsualEqRhs.hh>
0015 #include <Geant4/G4MagneticField.hh>
0016 #include <Geant4/G4PVPlacement.hh>
0017 #include <Geant4/G4PhysicalConstants.hh>
0018 #include <Geant4/G4QuadrupoleMagField.hh>
0019 #include <Geant4/G4RotationMatrix.hh>
0020 #include <Geant4/G4String.hh>  // for G4String
0021 #include <Geant4/G4SystemOfUnits.hh>
0022 #include <Geant4/G4ThreeVector.hh>  // for G4ThreeVector
0023 #include <Geant4/G4Transform3D.hh>  // for G4Transform3D
0024 #include <Geant4/G4Tubs.hh>
0025 #include <Geant4/G4Types.hh>  // for G4double, G4bool
0026 #include <Geant4/G4UniformMagField.hh>
0027 #include <Geant4/G4VisAttributes.hh>
0028
0029 #include <CLHEP/Units/SystemOfUnits.h>  // for cm, deg, tesla, twopi, meter
0030
0031 #include <cstdlib>   // for exit
0032 #include <iostream>  // for operator<<, basic_ostream
0033
0034 class G4Material;
0035 class G4VSolid;
0036 class PHCompositeNode;
0037 class PHG4Subsystem;
0038
0039 using namespace std;
0040
0041 //_______________________________________________________________
0042 PHG4BeamlineMagnetDetector::PHG4BeamlineMagnetDetector(PHG4Subsystem *subsys, PHCompositeNode *Node, PHParameters *parameters, const std::string &dnam, const int lyr)
0043   : PHG4Detector(subsys, Node, dnam)
0044   , params(parameters)
0045   , magnet_physi(nullptr)
0046   , cylinder_physi(nullptr)
0047   , layer(lyr)
0048 {
0049 }
0050
0051 //_______________________________________________________________
0052 bool PHG4BeamlineMagnetDetector::IsInBeamlineMagnet(const G4VPhysicalVolume *volume) const
0053 {
0054   if (volume == magnet_physi)
0055   {
0056     return true;
0057   }
0058   return false;
0059 }
0060
0061 //_______________________________________________________________
0062 void PHG4BeamlineMagnetDetector::ConstructMe(G4LogicalVolume *logicMother)
0063 {
0064   G4Material *TrackerMaterial = GetDetectorMaterial(params->get_string_param("material"));
0065
0066   G4VisAttributes *fieldVis = new G4VisAttributes();
0067   PHG4Utils::SetColour(fieldVis, "BlackHole");
0068   fieldVis->SetVisibility(false);
0069   fieldVis->SetForceSolid(false);
0070
0071   G4VisAttributes *siliconVis = new G4VisAttributes();
0072   if (params->get_int_param("blackhole"))
0073   {
0074     PHG4Utils::SetColour(siliconVis, "BlackHole");
0075     siliconVis->SetVisibility(false);
0076     siliconVis->SetForceSolid(false);
0077   }
0078   else
0079   {
0080     PHG4Utils::SetColour(siliconVis, params->get_string_param("material"));
0081     siliconVis->SetVisibility(true);
0082     siliconVis->SetForceSolid(true);
0083   }
0084
0085   /* Define origin vector (center of magnet) */
0086   G4ThreeVector origin(params->get_double_param("place_x") * cm,
0087                        params->get_double_param("place_y") * cm,
0088                        params->get_double_param("place_z") * cm);
0089
0090   /* Define magnet rotation matrix */
0091   G4RotationMatrix *rotm = new G4RotationMatrix();
0092   rotm->rotateX(params->get_double_param("rot_x") * deg);
0093   rotm->rotateY(params->get_double_param("rot_y") * deg);
0094   rotm->rotateZ(params->get_double_param("rot_z") * deg);
0095
0096   /* Creating a magnetic field */
0097   G4MagneticField *magField = nullptr;
0098
0099   string magnettype = params->get_string_param("magtype");
0100   if (magnettype == "dipole")
0101   {
0102     G4double fieldValue = params->get_double_param("field_y") * tesla;
0103     magField = new G4UniformMagField(G4ThreeVector(0., fieldValue, 0.));
0104
0105     if (Verbosity() > 0)
0106     {
0107       cout << "Creating DIPOLE with field " << fieldValue << " and name " << GetName() << endl;
0108     }
0109   }
0110   else if (magnettype == "quadrupole")
0111   {
0112     G4double fieldGradient = params->get_double_param("fieldgradient") * tesla / meter;
0113
0114     /* G4MagneticField::GetFieldValue( pos*, B* ) uses GLOBAL coordinates, not local.
0115      * Therefore, place magnetic field center at the correct location and angle for the
0116      * magnet AND do the same transformations for the logical volume (see below). */
0117     magField = new G4QuadrupoleMagField(fieldGradient, origin, rotm);
0118     //      magField = new PHG4QuadrupoleMagField ( fieldGradient, origin, rotm );
0119
0120     if (Verbosity() > 0)
0121     {
0122       cout << "Creating QUADRUPOLE with gradient " << fieldGradient << " and name " << GetName() << endl;
0123       cout << "at x, y, z = " << origin.x() << " , " << origin.y() << " , " << origin.z() << endl;
0124       cout << "with rotation around x, y, z axis  by: " << rotm->phiX() << ", " << rotm->phiY() << ", " << rotm->phiZ() << endl;
0125     }
0126   }
0127
0128   if (!magField)
0129   {
0130     cout << PHWHERE << " Aborting, no magnetic field specified for " << GetName() << endl;
0131     exit(1);
0132   }
0133
0134   /* Set up Geant4 field manager */
0135   G4Mag_UsualEqRhs *localEquation = new G4Mag_UsualEqRhs(magField);
0136   G4ClassicalRK4 *localStepper = new G4ClassicalRK4(localEquation);
0137   G4double minStep = 0.25 * mm;  // minimal step, 1 mm is default
0138   G4ChordFinder *localChordFinder = new G4ChordFinder(magField, minStep, localStepper);
0139
0140   G4FieldManager *fieldMgr = new G4FieldManager();
0141   fieldMgr->SetDetectorField(magField);
0142   fieldMgr->SetChordFinder(localChordFinder);
0143
0144   /* Add volume with magnetic field */
0145   double radius = params->get_double_param("radius") * cm;
0146   double thickness = params->get_double_param("thickness") * cm;
0147
0148   G4VSolid *magnet_solid = new G4Tubs(GetName(),
0149                                       0,
0150                                       radius + thickness,
0151                                       params->get_double_param("length") * cm / 2., 0, twopi);
0152
0153   G4LogicalVolume *magnet_logic = new G4LogicalVolume(magnet_solid,
0154                                                       GetDetectorMaterial("G4_Galactic"),
0155                                                       GetName(),
0156                                                       nullptr, nullptr, nullptr);
0157   magnet_logic->SetVisAttributes(fieldVis);
0158
0159   /* Set field manager for logical volume */
0160   G4bool allLocal = true;
0161   magnet_logic->SetFieldManager(fieldMgr, allLocal);
0162
0163   /* create magnet physical volume */
0164   magnet_physi = new G4PVPlacement(G4Transform3D(*rotm,
0165                                                  G4ThreeVector(params->get_double_param("place_x") * cm,
0166                                                                params->get_double_param("place_y") * cm,
0167                                                                params->get_double_param("place_z") * cm)),
0168                                    magnet_logic,
0169                                    GetName(),
0170                                    logicMother, false, false, OverlapCheck());
0171
0172   /* Add volume with solid magnet material */
0173   G4VSolid *cylinder_solid = new G4Tubs(G4String(GetName().append("_Solid")),
0174                                         radius,
0175                                         radius + thickness,
0176                                         params->get_double_param("length") * cm / 2., 0, twopi);
0177   G4LogicalVolume *cylinder_logic = new G4LogicalVolume(cylinder_solid,
0178                                                         TrackerMaterial,
0179                                                         G4String(GetName()),
0180                                                         nullptr, nullptr, nullptr);
0181   cylinder_logic->SetVisAttributes(siliconVis);
0182
0183   cylinder_physi = new G4PVPlacement(nullptr, G4ThreeVector(0, 0, 0),
0184                                      cylinder_logic,
0185                                      G4String(GetName().append("_Solid")),
0186                                      magnet_logic, false, false, OverlapCheck());
0187 }