// // ******************************************************************** // * License and Disclaimer * // * * // * The Geant4 software is copyright of the Copyright Holders of * // * the Geant4 Collaboration. It is provided under the terms and * // * conditions of the Geant4 Software License, included in the file * // * LICENSE and available at http://cern.ch/geant4/license . These * // * include a list of copyright holders. * // * * // * Neither the authors of this software system, nor their employing * // * institutes,nor the agencies providing financial support for this * // * work make any representation or warranty, express or implied, * // * regarding this software system or assume any liability for its * // * use. Please see the license in the file LICENSE and URL above * // * for the full disclaimer and the limitation of liability. * // * * // * This code implementation is the result of the scientific and * // * technical work of the GEANT4 collaboration. * // * By using, copying, modifying or distributing the software (or * // * any work based on the software) you agree to acknowledge its * // * use in resulting scientific publications, and indicate your * // * acceptance of all terms of the Geant4 Software license. * // ******************************************************************** // /// \file DetectorConstruction.cc /// \brief Implementation of the B5::DetectorConstruction class #include "DetectorConstruction.hh" #include "CellParameterisation.hh" #include "Constants.hh" #include "DriftChamberSD.hh" #include "EmCalorimeterSD.hh" #include "HadCalorimeterSD.hh" #include "HodoscopeSD.hh" #include "MagneticField.hh" #include "G4Box.hh" #include "G4Colour.hh" #include "G4FieldManager.hh" #include "G4GenericMessenger.hh" #include "G4LogicalVolume.hh" #include "G4Material.hh" #include "G4NistManager.hh" #include "G4PVParameterised.hh" #include "G4PVPlacement.hh" #include "G4PVReplica.hh" #include "G4RotationMatrix.hh" #include "G4RunManager.hh" #include "G4SDManager.hh" #include "G4SystemOfUnits.hh" #include "G4ThreeVector.hh" #include "G4Tubs.hh" #include "G4UserLimits.hh" #include "G4VisAttributes.hh" namespace B5 { //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... G4ThreadLocal MagneticField* DetectorConstruction::fMagneticField = nullptr; G4ThreadLocal G4FieldManager* DetectorConstruction::fFieldMgr = nullptr; //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... DetectorConstruction::DetectorConstruction() { fArmRotation = new G4RotationMatrix(); fArmRotation->rotateY(fArmAngle); // define commands for this class DefineCommands(); } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... DetectorConstruction::~DetectorConstruction() { delete fArmRotation; delete fMessenger; } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... G4VPhysicalVolume* DetectorConstruction::Construct() { // Construct materials ConstructMaterials(); auto air = G4Material::GetMaterial("G4_AIR"); // auto argonGas = G4Material::GetMaterial("_Ar"); auto argonGas = G4Material::GetMaterial("G4_Ar"); auto scintillator = G4Material::GetMaterial("G4_PLASTIC_SC_VINYLTOLUENE"); auto csI = G4Material::GetMaterial("G4_CESIUM_IODIDE"); auto lead = G4Material::GetMaterial("G4_Pb"); // Option to switch on/off checking of volumes overlaps // G4bool checkOverlaps = true; // geometries -------------------------------------------------------------- // experimental hall (world volume) auto worldSolid = new G4Box("worldBox", 10. * m, 3. * m, 10. * m); auto worldLogical = new G4LogicalVolume(worldSolid, air, "worldLogical"); auto worldPhysical = new G4PVPlacement(nullptr, G4ThreeVector(), worldLogical, "worldPhysical", nullptr, false, 0, checkOverlaps); // Tube with Local Magnetic field auto magneticSolid = new G4Tubs("magneticTubs", 0., 1. * m, 1. * m, 0., 360. * deg); fMagneticLogical = new G4LogicalVolume(magneticSolid, air, "magneticLogical"); // placement of Tube auto fieldRot = new G4RotationMatrix(); fieldRot->rotateX(90. * deg); new G4PVPlacement(fieldRot, G4ThreeVector(), fMagneticLogical, "magneticPhysical", worldLogical, false, 0, checkOverlaps); // set step limit in tube with magnetic field auto userLimits = new G4UserLimits(1 * m); fMagneticLogical->SetUserLimits(userLimits); // first arm auto firstArmSolid = new G4Box("firstArmBox", 1.5 * m, 1. * m, 3. * m); auto firstArmLogical = new G4LogicalVolume(firstArmSolid, air, "firstArmLogical"); new G4PVPlacement(nullptr, G4ThreeVector(0., 0., -5. * m), firstArmLogical, "firstArmPhysical", worldLogical, false, 0, checkOverlaps); // second arm auto secondArmSolid = new G4Box("secondArmBox", 2. * m, 2. * m, 3.5 * m); auto secondArmLogical = new G4LogicalVolume(secondArmSolid, air, "secondArmLogical"); auto x = -5. * m * std::sin(fArmAngle); auto z = 5. * m * std::cos(fArmAngle); fSecondArmPhys = new G4PVPlacement(fArmRotation, G4ThreeVector(x, 0., z), secondArmLogical, "fSecondArmPhys", worldLogical, false, 0, checkOverlaps); // hodoscopes in first arm auto hodoscope1Solid = new G4Box("hodoscope1Box", 5. * cm, 20. * cm, 0.5 * cm); fHodoscope1Logical = new G4LogicalVolume(hodoscope1Solid, scintillator, "hodoscope1Logical"); for (auto i = 0; i < kNofHodoscopes1; i++) { G4double x1 = (i - kNofHodoscopes1 / 2) * 10. * cm; new G4PVPlacement(nullptr, G4ThreeVector(x1, 0., -1.5 * m), fHodoscope1Logical, "hodoscope1Physical", firstArmLogical, false, i, checkOverlaps); } // drift chambers in first arm auto chamber1Solid = new G4Box("chamber1Box", 1. * m, 30. * cm, 1. * cm); auto chamber1Logical = new G4LogicalVolume(chamber1Solid, argonGas, "chamber1Logical"); for (auto i = 0; i < kNofChambers; i++) { G4double z1 = (i - kNofChambers / 2) * 0.5 * m; new G4PVPlacement(nullptr, G4ThreeVector(0., 0., z1), chamber1Logical, "chamber1Physical", firstArmLogical, false, i, checkOverlaps); } // "virtual" wire plane auto wirePlane1Solid = new G4Box("wirePlane1Box", 1. * m, 30. * cm, 0.1 * mm); fWirePlane1Logical = new G4LogicalVolume(wirePlane1Solid, argonGas, "wirePlane1Logical"); new G4PVPlacement(nullptr, G4ThreeVector(0., 0., 0.), fWirePlane1Logical, "wirePlane1Physical", chamber1Logical, false, 0, checkOverlaps); // hodoscopes in second arm auto hodoscope2Solid = new G4Box("hodoscope2Box", 5. * cm, 20. * cm, 0.5 * cm); fHodoscope2Logical = new G4LogicalVolume(hodoscope2Solid, scintillator, "hodoscope2Logical"); for (auto i = 0; i < kNofHodoscopes2; i++) { G4double x2 = (i - kNofHodoscopes2 / 2) * 10. * cm; new G4PVPlacement(nullptr, G4ThreeVector(x2, 0., 0.), fHodoscope2Logical, "hodoscope2Physical", secondArmLogical, false, i, checkOverlaps); } // drift chambers in second arm auto chamber2Solid = new G4Box("chamber2Box", 1.5 * m, 30. * cm, 1. * cm); auto chamber2Logical = new G4LogicalVolume(chamber2Solid, argonGas, "chamber2Logical"); for (auto i = 0; i < kNofChambers; i++) { G4double z2 = (i - kNofChambers / 2) * 0.5 * m - 1.5 * m; new G4PVPlacement(nullptr, G4ThreeVector(0., 0., z2), chamber2Logical, "chamber2Physical", secondArmLogical, false, i, checkOverlaps); } // "virtual" wire plane auto wirePlane2Solid = new G4Box("wirePlane2Box", 1.5 * m, 30. * cm, 0.1 * mm); fWirePlane2Logical = new G4LogicalVolume(wirePlane2Solid, argonGas, "wirePlane2Logical"); new G4PVPlacement(nullptr, G4ThreeVector(0., 0., 0.), fWirePlane2Logical, "wirePlane2Physical", chamber2Logical, false, 0, checkOverlaps); // CsI calorimeter auto emCalorimeterSolid = new G4Box("EMcalorimeterBox", 1.5 * m, 30. * cm, 15. * cm); auto emCalorimeterLogical = new G4LogicalVolume(emCalorimeterSolid, csI, "EMcalorimeterLogical"); new G4PVPlacement(nullptr, G4ThreeVector(0., 0., 2. * m), emCalorimeterLogical, "EMcalorimeterPhysical", secondArmLogical, false, 0, checkOverlaps); // EMcalorimeter cells auto cellSolid = new G4Box("cellBox", 7.5 * cm, 7.5 * cm, 15. * cm); fCellLogical = new G4LogicalVolume(cellSolid, csI, "cellLogical"); G4VPVParameterisation* cellParam = new CellParameterisation(); new G4PVParameterised("cellPhysical", fCellLogical, emCalorimeterLogical, kXAxis, kNofEmCells, cellParam); // hadron calorimeter auto hadCalorimeterSolid = new G4Box("HadCalorimeterBox", 1.5 * m, 30. * cm, 50. * cm); auto hadCalorimeterLogical = new G4LogicalVolume(hadCalorimeterSolid, lead, "HadCalorimeterLogical"); new G4PVPlacement(nullptr, G4ThreeVector(0., 0., 3. * m), hadCalorimeterLogical, "HadCalorimeterPhysical", secondArmLogical, false, 0, checkOverlaps); // hadron calorimeter column auto HadCalColumnSolid = new G4Box("HadCalColumnBox", 15. * cm, 30. * cm, 50. * cm); auto HadCalColumnLogical = new G4LogicalVolume(HadCalColumnSolid, lead, "HadCalColumnLogical"); new G4PVReplica("HadCalColumnPhysical", HadCalColumnLogical, hadCalorimeterLogical, kXAxis, kNofHadColumns, 30. * cm); // hadron calorimeter cell auto HadCalCellSolid = new G4Box("HadCalCellBox", 15. * cm, 15. * cm, 50. * cm); auto HadCalCellLogical = new G4LogicalVolume(HadCalCellSolid, lead, "HadCalCellLogical"); new G4PVReplica("HadCalCellPhysical", HadCalCellLogical, HadCalColumnLogical, kYAxis, kNofHadRows, 30. * cm); // hadron calorimeter layers auto HadCalLayerSolid = new G4Box("HadCalLayerBox", 15. * cm, 15. * cm, 2.5 * cm); auto HadCalLayerLogical = new G4LogicalVolume(HadCalLayerSolid, lead, "HadCalLayerLogical"); new G4PVReplica("HadCalLayerPhysical", HadCalLayerLogical, HadCalCellLogical, kZAxis, kNofHadCells, 5. * cm); // scintillator plates auto HadCalScintiSolid = new G4Box("HadCalScintiBox", 15. * cm, 15. * cm, 0.5 * cm); fHadCalScintiLogical = new G4LogicalVolume(HadCalScintiSolid, scintillator, "HadCalScintiLogical"); new G4PVPlacement(nullptr, G4ThreeVector(0., 0., 2. * cm), fHadCalScintiLogical, "HadCalScintiPhysical", HadCalLayerLogical, false, 0, checkOverlaps); // visualization attributes ------------------------------------------------ G4VisAttributes invisible(G4VisAttributes::GetInvisible()); G4VisAttributes invisibleBlue(false, G4Colour::Blue()); G4VisAttributes invisibleGreen(false, G4Colour::Green()); G4VisAttributes invisibleYellow(false, G4Colour::Yellow()); G4VisAttributes blue(G4Colour::Blue()); G4VisAttributes cgray(G4Colour::Gray()); G4VisAttributes green(G4Colour::Green()); G4VisAttributes red(G4Colour::Red()); G4VisAttributes yellow(G4Colour::Yellow()); worldLogical->SetVisAttributes(invisible); firstArmLogical->SetVisAttributes(invisible); secondArmLogical->SetVisAttributes(invisible); fMagneticLogical->SetVisAttributes(cgray); fHodoscope1Logical->SetVisAttributes(red); fHodoscope2Logical->SetVisAttributes(red); chamber1Logical->SetVisAttributes(green); chamber2Logical->SetVisAttributes(green); fWirePlane1Logical->SetVisAttributes(invisibleGreen); fWirePlane2Logical->SetVisAttributes(invisibleGreen); emCalorimeterLogical->SetVisAttributes(invisibleYellow); fCellLogical->SetVisAttributes(yellow); hadCalorimeterLogical->SetVisAttributes(blue); HadCalColumnLogical->SetVisAttributes(invisibleBlue); HadCalCellLogical->SetVisAttributes(invisibleBlue); HadCalLayerLogical->SetVisAttributes(invisibleBlue); fHadCalScintiLogical->SetVisAttributes(invisibleBlue); // return the world physical volume ---------------------------------------- return worldPhysical; } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void DetectorConstruction::ConstructSDandField() { // sensitive detectors ----------------------------------------------------- auto sdManager = G4SDManager::GetSDMpointer(); G4String SDname; auto hodoscope1 = new HodoscopeSD(SDname = "/hodoscope1"); sdManager->AddNewDetector(hodoscope1); fHodoscope1Logical->SetSensitiveDetector(hodoscope1); auto hodoscope2 = new HodoscopeSD(SDname = "/hodoscope2"); sdManager->AddNewDetector(hodoscope2); fHodoscope2Logical->SetSensitiveDetector(hodoscope2); auto chamber1 = new DriftChamberSD(SDname = "/chamber1"); sdManager->AddNewDetector(chamber1); fWirePlane1Logical->SetSensitiveDetector(chamber1); auto chamber2 = new DriftChamberSD(SDname = "/chamber2"); sdManager->AddNewDetector(chamber2); fWirePlane2Logical->SetSensitiveDetector(chamber2); auto emCalorimeter = new EmCalorimeterSD(SDname = "/EMcalorimeter"); sdManager->AddNewDetector(emCalorimeter); fCellLogical->SetSensitiveDetector(emCalorimeter); auto hadCalorimeter = new HadCalorimeterSD(SDname = "/HadCalorimeter"); sdManager->AddNewDetector(hadCalorimeter); fHadCalScintiLogical->SetSensitiveDetector(hadCalorimeter); // magnetic field ---------------------------------------------------------- fMagneticField = new MagneticField(); fFieldMgr = new G4FieldManager(); fFieldMgr->SetDetectorField(fMagneticField); fFieldMgr->CreateChordFinder(fMagneticField); G4bool forceToAllDaughters = true; fMagneticLogical->SetFieldManager(fFieldMgr, forceToAllDaughters); } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void DetectorConstruction::ConstructMaterials() { auto nistManager = G4NistManager::Instance(); // Air nistManager->FindOrBuildMaterial("G4_AIR"); // Argon gas nistManager->FindOrBuildMaterial("G4_Ar"); // With a density different from the one defined in NIST // G4double density = 1.782e-03*g/cm3; // nistManager->BuildMaterialWithNewDensity("_Ar","G4_Ar",density); // !! cases segmentation fault // Scintillator // (PolyVinylToluene, C_9H_10) nistManager->FindOrBuildMaterial("G4_PLASTIC_SC_VINYLTOLUENE"); // CsI nistManager->FindOrBuildMaterial("G4_CESIUM_IODIDE"); // Lead nistManager->FindOrBuildMaterial("G4_Pb"); // Vacuum "Galactic" // nistManager->FindOrBuildMaterial("G4_Galactic"); // Vacuum "Air with low density" // auto air = G4Material::GetMaterial("G4_AIR"); // G4double density = 1.0e-5*air->GetDensity(); // nistManager // ->BuildMaterialWithNewDensity("Air_lowDensity", "G4_AIR", density); G4cout << G4endl << "The materials defined are : " << G4endl << G4endl; G4cout << *(G4Material::GetMaterialTable()) << G4endl; } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void DetectorConstruction::SetArmAngle(G4double val) { if (!fSecondArmPhys) { G4cerr << "Detector has not yet been constructed." << G4endl; return; } fArmAngle = val; *fArmRotation = G4RotationMatrix(); // make it unit vector fArmRotation->rotateY(fArmAngle); auto x = -5. * m * std::sin(fArmAngle); auto z = 5. * m * std::cos(fArmAngle); fSecondArmPhys->SetTranslation(G4ThreeVector(x, 0., z)); // tell G4RunManager that we change the geometry G4RunManager::GetRunManager()->GeometryHasBeenModified(); } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void DetectorConstruction::DefineCommands() { // Define /B5/detector command directory using generic messenger class fMessenger = new G4GenericMessenger(this, "/B5/detector/", "Detector control"); // armAngle command auto& armAngleCmd = fMessenger->DeclareMethodWithUnit( "armAngle", "deg", &DetectorConstruction::SetArmAngle, "Set rotation angle of the second arm."); armAngleCmd.SetParameterName("angle", true); armAngleCmd.SetRange("angle>=0. && angle<180."); armAngleCmd.SetDefaultValue("30."); } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... } // namespace B5