GB03/src/GB03DetectorConstruction.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 /// \file GB03DetectorConstruction.cc
0028 /// \brief Implementation of the GB03DetectorConstruction class
0029
0030 #include "GB03DetectorConstruction.hh"
0031
0032 #include "GB03BOptrGeometryBasedBiasing.hh"
0033 #include "GB03DetectorMessenger.hh"
0034
0035 #include "G4Box.hh"
0036 #include "G4Colour.hh"
0037 #include "G4LogicalVolume.hh"
0038 #include "G4Material.hh"
0039 #include "G4MultiFunctionalDetector.hh"
0040 #include "G4PSEnergyDeposit.hh"
0041 #include "G4PSFlatSurfaceFlux.hh"
0042 #include "G4PVPlacement.hh"
0043 #include "G4PVReplica.hh"
0044 #include "G4PhysicalConstants.hh"
0045 #include "G4RunManager.hh"
0046 #include "G4SDChargedFilter.hh"
0047 #include "G4SDManager.hh"
0048 #include "G4SDNeutralFilter.hh"
0049 #include "G4SystemOfUnits.hh"
0050 #include "G4VisAttributes.hh"
0051 #include "G4ios.hh"
0052
0053 G4int GB03DetectorConstruction::fNumberOfLayers = 40;
0054 G4ThreadLocal G4bool GB03DetectorConstruction::fConstructedSDandField = false;
0055
0056 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0057
0058 GB03DetectorConstruction::GB03DetectorConstruction(G4bool bf)
0059   : G4VUserDetectorConstruction(),
0060     fTotalThickness(2.0 * m),
0061     fLayerThickness(0.),
0062     fConstructed(false),
0063     fWorldMaterial(nullptr),
0064     fAbsorberMaterial(nullptr),
0065     fGapMaterial(nullptr),
0066     fLayerSolid(nullptr),
0067     fGapSolid(nullptr),
0068     fWorldLogical(nullptr),
0069     fCalorLogical(nullptr),
0070     fLayerLogical(nullptr),
0071     fGapLogical(nullptr),
0072     fWorldPhysical(nullptr),
0073     fCalorPhysical(nullptr),
0074     fLayerPhysical(nullptr),
0075     fGapPhysical(nullptr),
0076     fDetectorMessenger(nullptr),
0077     fVerboseLevel(1),
0078     fBiasingFlag(bf)
0079 {
0080   fLayerThickness = fTotalThickness / fNumberOfLayers;
0081   fCalName = "Calor";
0082   fDetectorMessenger = new GB03DetectorMessenger(this);
0083 }
0084
0085 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0086
0087 GB03DetectorConstruction::~GB03DetectorConstruction()
0088 {
0089   delete fDetectorMessenger;
0090 }
0091
0092 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0093
0094 G4VPhysicalVolume* GB03DetectorConstruction::Construct()
0095 {
0096   if (!fConstructed) {
0097     fConstructed = true;
0098     DefineMaterials();
0099     SetupGeometry();
0100   }
0101   if (GetVerboseLevel() > 0) {
0102     PrintCalorParameters();
0103   }
0104
0105   return fWorldPhysical;
0106 }
0107
0108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0109
0110 void GB03DetectorConstruction::ConstructSDandField()
0111 {
0112   if (!fConstructedSDandField) {
0113     fConstructedSDandField = true;
0114     SetupDetectors();
0115     if (fBiasingFlag) {
0116       SetupBiasing();
0117     }
0118   }
0119 }
0120
0121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0122
0123 void GB03DetectorConstruction::DefineMaterials()
0124 {
0125   G4String name, symbol;  // a=mass of a mole;
0126   G4double a, z, density;  // z=mean number of protons;
0127   G4int iz;  // iz=number of protons  in an isotope;
0128   G4int n;  // n=number of nucleons in an isotope;
0129
0130   G4int ncomponents, natoms;
0131   G4double abundance, fractionmass;
0132   G4double temperature, pressure;
0133
0134   //
0135   // define Elements
0136   //
0137
0138   a = 1.01 * g / mole;
0139   auto H = new G4Element(name = "Hydrogen", symbol = "H", z = 1., a);
0140
0141   a = 12.01 * g / mole;
0142   auto C = new G4Element(name = "Carbon", symbol = "C", z = 6., a);
0143
0144   a = 14.01 * g / mole;
0145   auto N = new G4Element(name = "Nitrogen", symbol = "N", z = 7., a);
0146
0147   a = 16.00 * g / mole;
0148   auto O = new G4Element(name = "Oxygen", symbol = "O", z = 8., a);
0149
0150   //
0151   // define an Element from isotopes, by relative abundance
0152   //
0153
0154   auto U5 = new G4Isotope(name = "U235", iz = 92, n = 235, a = 235.01 * g / mole);
0155   auto U8 = new G4Isotope(name = "U238", iz = 92, n = 238, a = 238.03 * g / mole);
0156
0157   auto U = new G4Element(name = "enriched Uranium", symbol = "U", ncomponents = 2);
0158   U->AddIsotope(U5, abundance = 90. * perCent);
0159   U->AddIsotope(U8, abundance = 10. * perCent);
0160
0161   //
0162   // define simple materials
0163   //
0164
0165   new G4Material(name = "Aluminium", z = 13., a = 26.98 * g / mole, density = 2.700 * g / cm3);
0166   new G4Material(name = "Silicon", z = 14., a = 28.09 * g / mole, density = 2.33 * g / cm3);
0167   new G4Material(name = "Iron", z = 26., a = 55.85 * g / mole, density = 7.87 * g / cm3);
0168   new G4Material(name = "ArgonGas", z = 18., a = 39.95 * g / mole, density = 1.782 * mg / cm3);
0169   new G4Material(name = "He", z = 2., a = 4.0 * g / mole, density = 0.1786e-03 * g / cm3);
0170
0171   density = 1.390 * g / cm3;
0172   a = 39.95 * g / mole;
0173   new G4Material(name = "liquidArgon", z = 18., a, density);
0174
0175   density = 11.35 * g / cm3;
0176   a = 207.19 * g / mole;
0177   auto Pb = new G4Material(name = "Lead", z = 82., a, density);
0178
0179   //
0180   // define a material from elements.   case 1: chemical molecule
0181   //
0182
0183   density = 1.000 * g / cm3;
0184   auto H2O = new G4Material(name = "Water", density, ncomponents = 2);
0185   H2O->AddElement(H, natoms = 2);
0186   H2O->AddElement(O, natoms = 1);
0187
0188   density = 1.032 * g / cm3;
0189   auto Sci = new G4Material(name = "Scintillator", density, ncomponents = 2);
0190   Sci->AddElement(C, natoms = 9);
0191   Sci->AddElement(H, natoms = 10);
0192
0193   //
0194   // define a material from elements.   case 2: mixture by fractional mass
0195   //
0196
0197   density = 1.290 * mg / cm3;
0198   auto Air = new G4Material(name = "Air", density, ncomponents = 2);
0199   Air->AddElement(N, fractionmass = 0.7);
0200   Air->AddElement(O, fractionmass = 0.3);
0201
0202   //
0203   // examples of vacuum
0204   //
0205
0206   density = universe_mean_density;
0207   pressure = 3.e-18 * pascal;
0208   temperature = 2.73 * kelvin;
0209   auto Vacuum = new G4Material(name = "Galactic", z = 1., a = 1.01 * g / mole, density, kStateGas,
0210                                temperature, pressure);
0211
0212   if (GetVerboseLevel() > 1) {
0213     G4cout << *(G4Material::GetMaterialTable()) << G4endl;
0214   }
0215
0216   // default materials of the calorimeter
0217   fWorldMaterial = Vacuum;
0218   fAbsorberMaterial = Pb;
0219   fGapMaterial = Sci;
0220 }
0221
0222 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0223
0224 void GB03DetectorConstruction::SetupGeometry()
0225 {
0226   //
0227   // World
0228   //
0229   G4VSolid* worldSolid = new G4Box("World", 2. * m, 2. * m, fTotalThickness * 2.);
0230   fWorldLogical = new G4LogicalVolume(worldSolid, fWorldMaterial, "World");
0231   fWorldPhysical =
0232     new G4PVPlacement(nullptr, G4ThreeVector(), fWorldLogical, "World", nullptr, false, 0);
0233
0234   //
0235   // Calorimeter
0236   //
0237   G4VSolid* calorSolid = new G4Box("Calor", 0.5 * m, 0.5 * m, fTotalThickness / 2.);
0238   fCalorLogical = new G4LogicalVolume(calorSolid, fAbsorberMaterial, fCalName);
0239   fCalorPhysical = new G4PVPlacement(nullptr, G4ThreeVector(0., 0., 0.), fCalorLogical, fCalName,
0240                                      fWorldLogical, false, 0);
0241
0242   //
0243   // Layers --- as absorbers
0244   //
0245   fLayerSolid = new G4Box("Layer", 0.5 * m, 0.5 * m, fLayerThickness / 2.);
0246   fLayerLogical = new G4LogicalVolume(fLayerSolid, fAbsorberMaterial, fCalName + "_LayerLog");
0247   fLayerPhysical = new G4PVReplica(fCalName + "_Layer", fLayerLogical, fCalorLogical, kZAxis,
0248                                    fNumberOfLayers, fLayerThickness);
0249
0250   //
0251   // Gap
0252   //
0253   fGapSolid = new G4Box("Gap", 0.5 * m, 0.5 * m, fLayerThickness / 4.);
0254   fGapLogical = new G4LogicalVolume(fGapSolid, fGapMaterial, fCalName + "_Gap");
0255   fGapPhysical = new G4PVPlacement(nullptr, G4ThreeVector(0., 0., fLayerThickness / 4.),
0256                                    fGapLogical, fCalName + "_gap", fLayerLogical, false, 0);
0257
0258   //
0259   // Visualization attributes
0260   //
0261   fWorldLogical->SetVisAttributes(G4VisAttributes::GetInvisible());
0262   auto simpleBoxVisAtt = new G4VisAttributes(G4Colour(1.0, 1.0, 1.0));
0263   simpleBoxVisAtt->SetVisibility(true);
0264   fCalorLogical->SetVisAttributes(simpleBoxVisAtt);
0265   fLayerLogical->SetVisAttributes(simpleBoxVisAtt);
0266   fGapLogical->SetVisAttributes(simpleBoxVisAtt);
0267 }
0268
0269 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0270
0271 void GB03DetectorConstruction::SetupDetectors()
0272 {
0273   G4SDManager::GetSDMpointer()->SetVerboseLevel(1);
0274   G4String filterName;
0275
0276   auto neutralFilter = new G4SDNeutralFilter(filterName = "neutralFilter");
0277   auto chargedFilter = new G4SDChargedFilter(filterName = "chargedFilter");
0278
0279   for (G4int j = 0; j < 2; j++) {
0280     // Loop counter j = 0 : absorber
0281     //                = 1 : gap
0282     G4String detName = fCalName;
0283     if (j == 0) {
0284       detName += "_abs";
0285     }
0286     else {
0287       detName += "_gap";
0288     }
0289     auto det = new G4MultiFunctionalDetector(detName);
0290     G4SDManager::GetSDMpointer()->AddNewDetector(det);
0291     // The second argument in each primitive means the "level" of geometrical
0292     // hierarchy, the copy number of that level is used as the key of the
0293     // G4THitsMap.
0294     // For absorber (j = 0), the copy number of its own physical volume is used.
0295     // For gap (j = 1), the copy number of its mother physical volume is used,
0296     // since there is only one physical volume of gap is placed with respect
0297     // to its mother.
0298     G4VPrimitiveScorer* primitive;
0299     primitive = new G4PSEnergyDeposit("eDep", j);
0300     det->RegisterPrimitive(primitive);
0301     primitive = new G4PSFlatSurfaceFlux("nNeutral", 1, j);
0302     primitive->SetFilter(neutralFilter);
0303     det->RegisterPrimitive(primitive);
0304     primitive = new G4PSFlatSurfaceFlux("nCharged", 1, j);
0305     primitive->SetFilter(chargedFilter);
0306     det->RegisterPrimitive(primitive);
0307
0308     if (j == 0) {
0309       SetSensitiveDetector(fLayerLogical, det);
0310     }
0311     else {
0312       SetSensitiveDetector(fGapLogical, det);
0313     }
0314   }
0315   G4SDManager::GetSDMpointer()->SetVerboseLevel(0);
0316 }
0317
0318 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0319
0320 void GB03DetectorConstruction::SetupBiasing()
0321 {
0322   auto biasingOperator = new GB03BOptrGeometryBasedBiasing();
0323   biasingOperator->AttachTo(fLayerLogical);
0324 }
0325
0326 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0327
0328 void GB03DetectorConstruction::PrintCalorParameters() const
0329 {
0330   G4cout << "--------------------------------------------------------" << G4endl;
0331   G4cout << " Absorber is made of " << fAbsorberMaterial->GetName() << G4endl << " Gap is made of "
0332          << fGapMaterial->GetName() << G4endl
0333          << "--------------------------------------------------------" << G4endl;
0334 }
0335
0336 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0337
0338 void GB03DetectorConstruction::SetAbsorberMaterial(G4String materialChoice)
0339 {
0340   // search the material by its name
0341   G4Material* pttoMaterial = G4Material::GetMaterial(materialChoice);
0342   if (pttoMaterial) {
0343     fAbsorberMaterial = pttoMaterial;
0344     if (fConstructed) {
0345       fCalorLogical->SetMaterial(fAbsorberMaterial);
0346       fLayerLogical->SetMaterial(fAbsorberMaterial);
0347     }
0348     G4RunManager::GetRunManager()->GeometryHasBeenModified();
0349     if (GetVerboseLevel() > 1) {
0350       PrintCalorParameters();
0351     }
0352   }
0353   else {
0354     G4cerr << materialChoice << " is not defined. - Command is ignored." << G4endl;
0355   }
0356 }
0357
0358 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0359
0360 G4String GB03DetectorConstruction::GetAbsorberMaterial() const
0361 {
0362   return fAbsorberMaterial->GetName();
0363 }
0364
0365 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0366
0367 void GB03DetectorConstruction::SetGapMaterial(G4String materialChoice)
0368 {
0369   // search the material by its name
0370   G4Material* pttoMaterial = G4Material::GetMaterial(materialChoice);
0371   if (pttoMaterial) {
0372     fGapMaterial = pttoMaterial;
0373     if (fConstructed) {
0374       fGapLogical->SetMaterial(fGapMaterial);
0375     }
0376     G4RunManager::GetRunManager()->GeometryHasBeenModified();
0377     if (GetVerboseLevel() > 1) {
0378       PrintCalorParameters();
0379     }
0380   }
0381   else {
0382     G4cerr << materialChoice << " is not defined. - Command is ignored." << G4endl;
0383   }
0384 }
0385
0386 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0387
0388 G4String GB03DetectorConstruction::GetGapMaterial() const
0389 {
0390   return fGapMaterial->GetName();
0391 }
0392
0393 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0394
0395 void GB03DetectorConstruction::SetNumberOfLayers(G4int nl)
0396 {
0397   fNumberOfLayers = nl;
0398   fLayerThickness = fTotalThickness / fNumberOfLayers;
0399   if (!fConstructed) return;
0400
0401   fLayerSolid->SetZHalfLength(fLayerThickness / 2.);
0402   fGapSolid->SetZHalfLength(fLayerThickness / 4.);
0403
0404   fCalorLogical->RemoveDaughter(fLayerPhysical);
0405   delete fLayerPhysical;
0406   fLayerPhysical = new G4PVReplica(fCalName + "_Layer", fLayerLogical, fCalorLogical, kZAxis,
0407                                    fNumberOfLayers, fLayerThickness);
0408   fGapPhysical->SetTranslation(G4ThreeVector(0., 0., fLayerThickness / 4.));
0409
0410   G4RunManager::GetRunManager()->GeometryHasBeenModified();
0411 }
0412
0413 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......