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