radiobiology/src/VoxelizedSensitiveDetector.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 radiobiology/src/VoxelizedSensitiveDetector.cc
0028 /// \brief Implementation of the RadioBio::VoxelizedSensitiveDetector class
0029
0030 #include "VoxelizedSensitiveDetector.hh"
0031
0032 #include "G4Box.hh"
0033 #include "G4LogicalVolume.hh"
0034 #include "G4PVReplica.hh"
0035 #include "G4PhysicalVolumeStore.hh"
0036 #include "G4RunManager.hh"
0037 #include "G4SDManager.hh"
0038 #include "G4SystemOfUnits.hh"
0039 #include "G4TransportationManager.hh"
0040 #include "G4VPhysicalVolume.hh"
0041
0042 #include "DetectorConstruction.hh"
0043 #include "SD.hh"
0044 #include "VoxelizedSensitiveDetectorMessenger.hh"
0045
0046 namespace RadioBio
0047 {
0048
0049 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0050
0051 VoxelizedSensitiveDetector* VoxelizedSensitiveDetector::fInstance = nullptr;
0052
0053 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0054
0055 VoxelizedSensitiveDetector* VoxelizedSensitiveDetector::CreateInstance(DetectorConstruction* det,
0056                                                                        double xWidth, double yWidth,
0057                                                                        double zWidth)
0058 {
0059   if (fInstance) {
0060     delete fInstance;
0061     G4Exception("VoxelizedSensitiveDetector::createInstance", "RecreatingVoxelization",
0062                 FatalException, "Creating another, new, instance of VoxelizedSensitiveDetector");
0063   }
0064   fInstance = new VoxelizedSensitiveDetector(det, xWidth, yWidth, zWidth);
0065   return fInstance;
0066 }
0067
0068 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0069
0070 VoxelizedSensitiveDetector* VoxelizedSensitiveDetector::GetInstance()
0071 {
0072   return fInstance;
0073 }
0074
0075 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0076
0077 VoxelizedSensitiveDetector::VoxelizedSensitiveDetector(DetectorConstruction* det, double xWidth,
0078                                                        double yWidth, double zWidth)
0079   : fDetector(det), fVoxelWidthX(xWidth), fVoxelWidthY(yWidth), fVoxelWidthZ(zWidth)
0080 {
0081   fVoxelizedSensitiveDetectorMessenger = new VoxelizedSensitiveDetectorMessenger(this);
0082   UpdateVoxelVolume();
0083   CalculateVoxelNumber();
0084 }
0085
0086 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0087
0088 VoxelizedSensitiveDetector::~VoxelizedSensitiveDetector()
0089 {
0090   delete fVoxelizedSensitiveDetectorMessenger;
0091 }
0092
0093 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0094
0095 void VoxelizedSensitiveDetector::UpdateVoxelVolume()
0096 {
0097   fVoxelVolume = fVoxelWidthX * fVoxelWidthY * fVoxelWidthZ;
0098   fVoxelDensity = fDetector->GetMaterial()->GetDensity();
0099   fVoxelMass = fVoxelVolume * fVoxelDensity;
0100 }
0101
0102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0103
0104 void VoxelizedSensitiveDetector::SetVoxelWidth(G4ThreeVector voxWidth)
0105 {
0106   fVoxelWidthX = voxWidth.getX();
0107   fVoxelWidthY = voxWidth.getY();
0108   fVoxelWidthZ = voxWidth.getZ();
0109   CalculateVoxelNumber();
0110 }
0111
0112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0113
0114 void VoxelizedSensitiveDetector::SetVoxelWidthX(G4double voxWidthX)
0115 {
0116   if (fVoxelWidthX == voxWidthX) return;
0117   fVoxelWidthX = voxWidthX;
0118   CalculateVoxelNumber();
0119 }
0120
0121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0122
0123 void VoxelizedSensitiveDetector::SetVoxelWidthY(G4double voxWidthY)
0124 {
0125   if (fVoxelWidthY == voxWidthY) return;
0126   fVoxelWidthY = voxWidthY;
0127   CalculateVoxelNumber();
0128 }
0129
0130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0131
0132 void VoxelizedSensitiveDetector::SetVoxelWidthZ(G4double voxWidthZ)
0133 {
0134   if (fVoxelWidthZ == voxWidthZ) return;
0135   fVoxelWidthZ = voxWidthZ;
0136   CalculateVoxelNumber();
0137 }
0138
0139 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0140
0141 // Calculte number of voxel approximating for an integer number of voxels
0142 // Then recalculates voxels size according to approximations
0143 void VoxelizedSensitiveDetector::CalculateVoxelNumber()
0144 {
0145   fVoxelNumberAlongX = G4int(fDetector->GetSizeX() / fVoxelWidthX);
0146   fVoxelWidthX = fDetector->GetSizeX() / G4double(fVoxelNumberAlongX);
0147
0148   fVoxelNumberAlongY = G4int(fDetector->GetSizeY() / fVoxelWidthY);
0149   fVoxelWidthY = fDetector->GetSizeY() / G4double(fVoxelNumberAlongY);
0150
0151   fVoxelNumberAlongZ = G4int(fDetector->GetSizeZ() / fVoxelWidthZ);
0152   fVoxelWidthZ = fDetector->GetSizeZ() / G4double(fVoxelNumberAlongZ);
0153
0154   if (fVoxelNumberAlongY % 2 == 0)
0155     G4Exception("VoxelizedSensitiveDetector::CalculateVoxelNumber", "VoxelNumberYEven", JustWarning,
0156                 "Trying to voxelize with an even number of voxels along the Y axis."
0157                 "Please select an odd number to prevent from warnings due to tracking");
0158
0159   if (fVoxelNumberAlongZ % 2 == 0)
0160     G4Exception("VoxelizedSensitiveDetector::CalculateVoxelNumber", "VoxelNumberZEven", JustWarning,
0161                 "Trying to voxelize with an even number of voxels along the Z axis."
0162                 "Please select an odd number to prevent from warnings due to tracking");
0163
0164   fTotalVoxelNumber = fVoxelNumberAlongX * fVoxelNumberAlongY * fVoxelNumberAlongZ;
0165
0166   UpdateVoxelVolume();
0167 }
0168
0169 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0170
0171 void VoxelizedSensitiveDetector::ConstructXDivision()
0172 {
0173   if (fWorldLogical == nullptr)
0174     G4Exception("VoxelizedSensitiveDetector::ConstructXDivision", "WorldNotInit", FatalException,
0175                 "Voxelizing without having a pointer to world logical volume!");
0176
0177   if (!fDetector)
0178     G4Exception("VoxelizedSensitiveDetector::ConstructXDivision", "DetConstInit", FatalException,
0179                 "Voxelizing without having a pointer to DetectorConstruction!");
0180
0181   fVoxelizedDetectorXDivision = new G4Box("VoxelizedDetectorXDivision", fVoxelWidthX / 2,
0182                                           fDetector->GetSizeY() / 2, fDetector->GetSizeZ() / 2);
0183
0184   fVoxelizedDetectorXDivisionLog =
0185     new G4LogicalVolume(fVoxelizedDetectorXDivision, fWorldLogical->GetMaterial(),
0186                         "VoxelizedDetectorXDivisionLog", 0, 0, 0);
0187
0188   fVoxelizedDetectorXDivisionPhys =
0189     new G4PVReplica("VoxelizedDetectorXDivisionPhys", fVoxelizedDetectorXDivisionLog, fWorldLogical,
0190                     kXAxis, fVoxelNumberAlongX, fVoxelWidthX);
0191 }
0192
0193 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0194
0195 void VoxelizedSensitiveDetector::ConstructYDivision()
0196 {
0197   fVoxelizedDetectorYDivision = new G4Box("VoxelizedDetectorYDivision", fVoxelWidthX / 2,
0198                                           fVoxelWidthY / 2, fDetector->GetSizeZ() / 2);
0199
0200   fVoxelizedDetectorYDivisionLog =
0201     new G4LogicalVolume(fVoxelizedDetectorYDivision, fWorldLogical->GetMaterial(),
0202                         "VoxelizedDetectorYDivisionLog", 0, 0, 0);
0203
0204   fVoxelizedDetectorYDivisionPhys =
0205     new G4PVReplica("VoxelizedDetectorYDivisionPhys", fVoxelizedDetectorYDivisionLog,
0206                     fVoxelizedDetectorXDivisionLog, kYAxis, fVoxelNumberAlongY, fVoxelWidthY);
0207 }
0208
0209 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0210
0211 void VoxelizedSensitiveDetector::ConstructZDivision()
0212 {
0213   fVoxelizedDetectorZDivision =
0214     new G4Box("VoxelizedDetectorZDivision", fVoxelWidthX / 2, fVoxelWidthY / 2, fVoxelWidthZ / 2);
0215
0216   fVoxelizedDetectorZDivisionLog =
0217     new G4LogicalVolume(fVoxelizedDetectorZDivision, fWorldLogical->GetMaterial(),
0218                         "VoxelizedDetectorZDivisionLog", 0, 0, 0);
0219
0220   fVoxelizedDetectorZDivisionPhys =
0221     new G4PVReplica("VoxelizedDetectorZDivisionPhys", fVoxelizedDetectorZDivisionLog,
0222                     fVoxelizedDetectorYDivisionPhys, kZAxis, fVoxelNumberAlongZ, fVoxelWidthZ);
0223
0224   fSensitiveLogicalVolume = fVoxelizedDetectorZDivisionLog;
0225 }
0226
0227 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0228
0229 // First voxelize along X, then Y, then Z
0230 G4bool VoxelizedSensitiveDetector::ConstructVoxelizedDetector()
0231 {
0232   // Creating X division
0233   ConstructXDivision();
0234
0235   // Creating Y division
0236   ConstructYDivision();
0237
0238   // Creating Z division
0239   ConstructZDivision();
0240
0241   // Set last, smallest volumes as sensitive
0242   fSensitiveLogicalVolume = fVoxelizedDetectorZDivisionLog;
0243   fIsBuilt = true;
0244
0245   return true;
0246 }
0247
0248 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0249
0250 void VoxelizedSensitiveDetector::UpdateVoxelizedGeometry()
0251 {
0252   // Nothing happens if the voxelized geometry is not built. But parameters are properly set.
0253   if (!fIsBuilt) {
0254     return;
0255   }
0256
0257   CalculateVoxelNumber();
0258
0259   // Volume that will be deleted in order to update
0260   G4VPhysicalVolume* myVol;
0261
0262   G4PhysicalVolumeStore* store = G4PhysicalVolumeStore::GetInstance();
0263
0264   myVol = store->GetVolume("VoxelizedDetectorXDivisionPhys");
0265   store->DeRegister(myVol);
0266   myVol = store->GetVolume("VoxelizedDetectorYDivisionPhys");
0267   store->DeRegister(myVol);
0268   myVol = store->GetVolume("VoxelizedDetectorZDivisionPhys");
0269   store->DeRegister(myVol);
0270   fVoxelizedDetectorXDivisionPhys =
0271     new G4PVReplica("VoxelizedDetectorXDivisionPhys", fVoxelizedDetectorXDivisionLog, fWorldLogical,
0272                     kXAxis, fVoxelNumberAlongX, fVoxelWidthX);
0273
0274   fVoxelizedDetectorYDivisionPhys =
0275     new G4PVReplica("VoxelizedDetectorYDivisionPhys", fVoxelizedDetectorYDivisionLog,
0276                     fVoxelizedDetectorXDivisionPhys, kYAxis, fVoxelNumberAlongY, fVoxelWidthY);
0277
0278   fVoxelizedDetectorZDivisionPhys =
0279     new G4PVReplica("VoxelizedDetectorZDivisionPhys", fVoxelizedDetectorZDivisionLog,
0280                     fVoxelizedDetectorYDivisionPhys, kZAxis, fVoxelNumberAlongZ, fVoxelWidthZ);
0281
0282   G4RunManager::GetRunManager()->GeometryHasBeenModified();
0283   G4RunManager::GetRunManager()->PhysicsHasBeenModified();
0284 }
0285
0286 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0287
0288 void VoxelizedSensitiveDetector::ConstructSD()
0289 {
0290   G4String sensitiveDetectorName = "VoxelizedDetector";
0291   G4String HCname = "LETdata";
0292
0293   SD* detectorSD = new SD(sensitiveDetectorName, HCname);
0294   G4SDManager::GetSDMpointer()->AddNewDetector(detectorSD);
0295   fSensitiveLogicalVolume->SetSensitiveDetector(detectorSD);
0296 }
0297
0298 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0299
0300 void VoxelizedSensitiveDetector::Construct()
0301 {
0302   ConstructVoxelizedDetector();
0303 }
0304
0305 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0306
0307 void VoxelizedSensitiveDetector::InitializeWorldPtr(G4VPhysicalVolume* pWorld)
0308 {
0309   if (pWorld == nullptr)
0310     G4Exception("VoxelizedSensitiveDetector::InitializeWorldPtr", "WorldinitNull", FatalException,
0311                 "Initializing Voxelization Class with a Null Pointer to World!");
0312   fWorldLogical = pWorld->GetLogicalVolume();
0313 }
0314
0315 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0316
0317 // Returns absolute voxel index given matrix indexes
0318 G4int VoxelizedSensitiveDetector::GetThisVoxelNumber(G4int x, G4int y, G4int z) const
0319 {
0320   G4int nz = GetVoxelNumberAlongZ();
0321   G4int ny = GetVoxelNumberAlongY();
0322
0323   return z + nz * (y + ny * (x));
0324 }
0325
0326 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0327
0328 }  // namespace RadioBio