TestEm7/src/DetectorConstruction.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 /// \file electromagnetic/TestEm7/src/DetectorConstruction.cc
0027 /// \brief Implementation of the DetectorConstruction class
0028 //
0029 //
0030 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0031 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0032
0033 #include "DetectorConstruction.hh"
0034
0035 #include "DetectorMessenger.hh"
0036
0037 #include "G4Box.hh"
0038 #include "G4FieldManager.hh"
0039 #include "G4GeometryManager.hh"
0040 #include "G4LogicalVolume.hh"
0041 #include "G4LogicalVolumeStore.hh"
0042 #include "G4Material.hh"
0043 #include "G4NistManager.hh"
0044 #include "G4PVPlacement.hh"
0045 #include "G4PhysicalConstants.hh"
0046 #include "G4PhysicalVolumeStore.hh"
0047 #include "G4RunManager.hh"
0048 #include "G4SolidStore.hh"
0049 #include "G4SystemOfUnits.hh"
0050 #include "G4TransportationManager.hh"
0051 #include "G4UniformMagField.hh"
0052 #include "G4UnitsTable.hh"
0053
0054 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0055
0056 DetectorConstruction::DetectorConstruction()
0057   : G4VUserDetectorConstruction(), fMagField(nullptr), fLAbsor(nullptr), fLWorld(nullptr)
0058 {
0059   // default parameter values
0060   fAbsorSizeX = fAbsorSizeYZ = 20 * cm;
0061   fWorldSizeX = fWorldSizeYZ = 1.2 * fAbsorSizeX;
0062
0063   fTallyNumber = 0;
0064   for (G4int j = 0; j < kMaxTally; j++) {
0065     fTallySize[j] = fTallyPosition[j] = G4ThreeVector(0., 0., 0.);
0066     fTallyMass[j] = 0.;
0067     fLTally[j] = nullptr;
0068   }
0069
0070   DefineMaterials();
0071
0072   // create commands for interactive definition of the detector
0073   fDetectorMessenger = new DetectorMessenger(this);
0074 }
0075
0076 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0077
0078 DetectorConstruction::~DetectorConstruction()
0079 {
0080   delete fDetectorMessenger;
0081 }
0082
0083 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0084
0085 void DetectorConstruction::DefineMaterials()
0086 {
0087   //
0088   // define Elements
0089   //
0090   G4double z, a;
0091
0092   G4Element* H = new G4Element("Hydrogen", "H", z = 1, a = 1.008 * g / mole);
0093   G4Element* N = new G4Element("Nitrogen", "N", z = 7, a = 14.01 * g / mole);
0094   G4Element* O = new G4Element("Oxygen", "O", z = 8, a = 16.00 * g / mole);
0095
0096   //
0097   // define Materials.
0098   //
0099   G4double density, temperature, pressure;
0100   G4int ncomponents, natoms;
0101   G4double fractionmass;
0102
0103   G4Material* H2O = new G4Material("Water", density = 1.0 * g / cm3, ncomponents = 2);
0104   H2O->AddElement(H, natoms = 2);
0105   H2O->AddElement(O, natoms = 1);
0106   H2O->GetIonisation()->SetMeanExcitationEnergy(78.0 * eV);
0107
0108   // In this line both G4_WATER and Water_1.05 will be constructed
0109   G4NistManager::Instance()->BuildMaterialWithNewDensity("Water_1.05", "G4_WATER", 1.05 * g / cm3);
0110
0111   G4Material* Air = new G4Material("Air", density = 1.290 * mg / cm3, ncomponents = 2);
0112   Air->AddElement(N, fractionmass = 0.7);
0113   Air->AddElement(O, fractionmass = 0.3);
0114
0115   density = 1.e-5 * g / cm3;
0116   pressure = 2.e-2 * bar;
0117   temperature = STP_Temperature;  // From PhysicalConstants.h .
0118   G4Material* vac = new G4Material("TechVacuum", density, 1, kStateGas, temperature, pressure);
0119   vac->AddMaterial(Air, 1.);
0120
0121   density = universe_mean_density;  // from PhysicalConstants.h
0122   pressure = 3.e-18 * pascal;
0123   temperature = 2.73 * kelvin;
0124   G4Material* vacuum = new G4Material("Galactic", z = 1, a = 1.008 * g / mole, density, kStateGas,
0125                                       temperature, pressure);
0126
0127   // default materials
0128   fAbsorMaterial = H2O;
0129   fWorldMaterial = vacuum;
0130 }
0131
0132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0133
0134 G4VPhysicalVolume* DetectorConstruction::Construct()
0135 {
0136   // World
0137   //
0138   G4Box* sWorld = new G4Box("World",  // name
0139                             fWorldSizeX / 2, fWorldSizeYZ / 2, fWorldSizeYZ / 2);  // dimensions
0140
0141   fLWorld = new G4LogicalVolume(sWorld,  // shape
0142                                 fWorldMaterial,  // material
0143                                 "World");  // name
0144
0145   G4VPhysicalVolume* pWorld = new G4PVPlacement(0,  // no rotation
0146                                                 G4ThreeVector(0., 0., 0.),  // at (0,0,0)
0147                                                 fLWorld,  // logical volume
0148                                                 "World",  // name
0149                                                 0,  // mother  volume
0150                                                 false,  // no boolean operation
0151                                                 0);  // copy number
0152   //
0153   // Absorber
0154   //
0155   G4Box* sAbsor = new G4Box("Absorber",  // name
0156                             fAbsorSizeX / 2, fAbsorSizeYZ / 2, fAbsorSizeYZ / 2);  // dimensions
0157
0158   fLAbsor = new G4LogicalVolume(sAbsor,  // shape
0159                                 fAbsorMaterial,  // material
0160                                 "Absorber");  // name
0161
0162   new G4PVPlacement(0,  // no rotation
0163                     G4ThreeVector(0., 0., 0.),  // at (0,0,0)
0164                     fLAbsor,  // logical volume
0165                     "Absorber",  // name
0166                     fLWorld,  // mother  volume
0167                     false,  // no boolean operation
0168                     0);  // copy number
0169   //
0170   // Tallies (optional)
0171   //
0172   if (fTallyNumber > 0) {
0173     for (G4int j = 0; j < fTallyNumber; ++j) {
0174       G4Box* sTally =
0175         new G4Box("Tally", fTallySize[j].x() / 2, fTallySize[j].y() / 2, fTallySize[j].z() / 2);
0176
0177       fLTally[j] = new G4LogicalVolume(sTally, fAbsorMaterial, "Tally");
0178
0179       new G4PVPlacement(0,  // no rotation
0180                         fTallyPosition[j],  // position
0181                         fLTally[j],  // logical volume
0182                         "Tally",  // name
0183                         fLAbsor,  // mother  volume
0184                         false,  // no boolean operation
0185                         j + 1);  // copy number
0186
0187       fTallyMass[j] =
0188         fTallySize[j].x() * fTallySize[j].y() * fTallySize[j].z() * (fAbsorMaterial->GetDensity());
0189     }
0190   }
0191
0192   PrintParameters();
0193
0194   //
0195   // always return the World volume
0196   //
0197   return pWorld;
0198 }
0199
0200 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0201
0202 void DetectorConstruction::PrintParameters() const
0203 {
0204   G4cout << *(G4Material::GetMaterialTable()) << G4endl;
0205   G4cout << "\n---------------------------------------------------------\n";
0206   G4cout << "---> The Absorber is " << G4BestUnit(fAbsorSizeX, "Length") << " of "
0207          << fAbsorMaterial->GetName() << G4endl;
0208   G4cout << "\n---------------------------------------------------------\n";
0209
0210   if (fTallyNumber > 0) {
0211     G4cout << "---> There are " << fTallyNumber << " tallies : " << G4endl;
0212     for (G4int j = 0; j < fTallyNumber; ++j) {
0213       G4cout << "fTally " << j << ": " << fAbsorMaterial->GetName()
0214              << ",  mass = " << G4BestUnit(fTallyMass[j], "Mass")
0215              << " size = " << G4BestUnit(fTallySize[j], "Length")
0216              << " position = " << G4BestUnit(fTallyPosition[j], "Length") << G4endl;
0217     }
0218     G4cout << "\n---------------------------------------------------------\n";
0219   }
0220 }
0221
0222 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0223
0224 void DetectorConstruction::SetSizeX(G4double value)
0225 {
0226   fAbsorSizeX = value;
0227   fWorldSizeX = 1.2 * fAbsorSizeX;
0228 }
0229
0230 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0231
0232 void DetectorConstruction::SetSizeYZ(G4double value)
0233 {
0234   fAbsorSizeYZ = value;
0235   fWorldSizeYZ = 1.2 * fAbsorSizeYZ;
0236 }
0237
0238 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0239
0240 void DetectorConstruction::SetMaterial(const G4String& materialChoice)
0241 {
0242   // search the material by its name
0243   G4Material* pttoMaterial = G4NistManager::Instance()->FindOrBuildMaterial(materialChoice);
0244   if (pttoMaterial && pttoMaterial != fAbsorMaterial) {
0245     // change target material everywhere
0246     fAbsorMaterial = pttoMaterial;
0247     for (G4int j = 0; j < fTallyNumber; ++j) {
0248       if (fLTally[j]) {
0249         fLTally[j]->SetMaterial(pttoMaterial);
0250         fTallyMass[j] =
0251           fTallySize[j].x() * fTallySize[j].y() * fTallySize[j].z() * (pttoMaterial->GetDensity());
0252       }
0253     }
0254     if (fLAbsor) {
0255       fLAbsor->SetMaterial(fAbsorMaterial);
0256       G4RunManager::GetRunManager()->PhysicsHasBeenModified();
0257     }
0258   }
0259 }
0260
0261 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0262
0263 void DetectorConstruction::SetWorldMaterial(const G4String& materialChoice)
0264 {
0265   // search the material by its name
0266   G4Material* pttoMaterial = G4NistManager::Instance()->FindOrBuildMaterial(materialChoice);
0267   if (pttoMaterial && pttoMaterial != fWorldMaterial) {
0268     fWorldMaterial = pttoMaterial;
0269     if (fLWorld) {
0270       fLWorld->SetMaterial(fAbsorMaterial);
0271       G4RunManager::GetRunManager()->PhysicsHasBeenModified();
0272     }
0273   }
0274 }
0275
0276 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0277
0278 void DetectorConstruction::SetMagField(G4double fieldValue)
0279 {
0280   // apply a global uniform magnetic field along Z axis
0281   G4FieldManager* fieldMgr = G4TransportationManager::GetTransportationManager()->GetFieldManager();
0282
0283   if (fMagField) delete fMagField;  // delete the existing magn field
0284
0285   if (fieldValue != 0.)  // create a new one if non nul
0286   {
0287     fMagField = new G4UniformMagField(G4ThreeVector(0., 0., fieldValue));
0288     fieldMgr->SetDetectorField(fMagField);
0289     fieldMgr->CreateChordFinder(fMagField);
0290   }
0291   else {
0292     fMagField = nullptr;
0293     fieldMgr->SetDetectorField(fMagField);
0294   }
0295 }
0296 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0297
0298 void DetectorConstruction::SetTallyNumber(G4int value)
0299 {
0300   if (value >= 0 && value < kMaxTally) {
0301     fTallyNumber = value;
0302   }
0303   else {
0304     G4cout << "### DetectorConstruction::SetTallyNumber WARNING: wrong tally "
0305            << "number " << value << " is ignored" << G4endl;
0306   }
0307 }
0308
0309 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0310
0311 void DetectorConstruction::SetTallySize(G4int j, const G4ThreeVector& value)
0312 {
0313   if (j >= 0 && j < kMaxTally) {
0314     fTallySize[j] = value;
0315   }
0316   else {
0317     G4cout << "### DetectorConstruction::SetTallyNumber WARNING: wrong tally "
0318            << "number " << j << " is ignored" << G4endl;
0319   }
0320 }
0321
0322 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0323
0324 void DetectorConstruction::SetTallyPosition(G4int j, const G4ThreeVector& value)
0325 {
0326   if (j >= 0 && j < kMaxTally) {
0327     fTallyPosition[j] = value;
0328   }
0329   else {
0330     G4cout << "### DetectorConstruction::SetTallyPosition WARNING: wrong tally "
0331            << "number " << j << " is ignored" << G4endl;
0332   }
0333 }
0334
0335 G4double DetectorConstruction::GetTallyMass(G4int j) const
0336 {
0337   if (j >= 0 && j < kMaxTally) {
0338     return fTallyMass[j];
0339   }
0340   else {
0341     G4cout << "### DetectorConstruction::GetTallyMass WARNING: wrong tally "
0342            << "number " << j << " is ignored" << G4endl;
0343     return 0.0;
0344   }
0345 }
0346
0347 const G4LogicalVolume* DetectorConstruction::GetLogicalTally(G4int j) const
0348 {
0349   if (j >= 0 && j < kMaxTally) {
0350     return fLTally[j];
0351   }
0352   else {
0353     G4cout << "### DetectorConstruction::GetLOgicalTally WARNING: wrong tally "
0354            << "number " << j << " is ignored" << G4endl;
0355     return nullptr;
0356   }
0357 }
0358
0359 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......