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