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