cellularPhantom/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 // -----------------------------------------------------------------------------
0027 //       MONTE CARLO SIMULATION OF REALISTIC GEOMETRY FROM MICROSCOPES IMAGES
0028 //
0029 // Authors and contributors:
0030 // P. Barberet (a), S. Incerti (a), N. H. Tran (a), L. Morelli (a,b)
0031 //
0032 // a) University of Bordeaux, CNRS, LP2i, UMR5797, Gradignan, France
0033 // b) Politecnico di Milano, Italy
0034 //
0035 // If you use this code, please cite the following publication:
0036 // P. Barberet et al.,
0037 // "Monte-Carlo dosimetry on a realistic cell monolayer
0038 // geometry exposed to alpha particles."
0039 // Ph. Barberet et al 2012 Phys. Med. Biol. 57 2189
0040 // doi: 110.1088/0031-9155/57/8/2189
0041 // -----------------------------------------------------------------------------
0042
0043 #include "DetectorConstruction.hh"
0044 #include "DetectorMessenger.hh"
0045
0046 #include "G4PhysicalConstants.hh"
0047 #include "G4NistManager.hh"
0048 #include "G4ProductionCuts.hh"
0049
0050 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0051
0052 DetectorConstruction::DetectorConstruction()
0053 :G4VUserDetectorConstruction()
0054 {
0055   fDetectorMessenger = new DetectorMessenger(this);
0056 }
0057
0058 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0059
0060 G4VPhysicalVolume *DetectorConstruction::Construct()
0061 {
0062   DefineMaterials();
0063   return ConstructLine();
0064 }
0065
0066 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0067
0068 void DetectorConstruction::DefineMaterials()
0069 {
0070   G4String name, symbol;
0071
0072   // Water and air are defined from NIST material database
0073   G4NistManager *man = G4NistManager::Instance();
0074
0075   G4Material *H2O = man->FindOrBuildMaterial("G4_WATER");
0076   G4Material *Air = man->FindOrBuildMaterial("G4_AIR");
0077
0078   fDefaultMaterial = Air;
0079   fPhantomMaterial = H2O;  // material is not relevant
0080                            // it will be changed by the ComputeMaterial
0081                            // method of the CellParameterisation
0082
0083   // Default materials
0084   if (fMediumMaterial == nullptr) {fMediumMaterial = H2O;}
0085   if (fRedMaterial == nullptr) {fRedMaterial = H2O;}
0086   if (fGreenMaterial == nullptr) {fGreenMaterial = H2O;}
0087   if (fBlueMaterial == nullptr) {fBlueMaterial = H2O;}
0088 }
0089
0090 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0091
0092 G4VPhysicalVolume *DetectorConstruction::ConstructLine() {
0093
0094   //*************
0095   // World volume
0096   //*************
0097
0098   fSolidWorld = new G4Box("World",                         //its name
0099                           fWorldSizeXY / 2, fWorldSizeXY / 2, fWorldSizeZ / 2);  //its size
0100
0101   fLogicWorld = new G4LogicalVolume(fSolidWorld,           //its solid
0102                                     fDefaultMaterial,      //its material
0103                                     "World");              //its name
0104
0105   fPhysiWorld = new G4PVPlacement(nullptr,                 //no rotation
0106                                   G4ThreeVector(),         //at (0,0,0)
0107                                   "World",                 //its name
0108                                   fLogicWorld,             //its logical volume
0109                                   nullptr,                 //its mother  volume
0110                                   false,                   //no boolean operation
0111                                   0);                      //copy number
0112
0113   //********************
0114   // Cell culture medium
0115   //********************
0116
0117   fSolidMedium = new G4Box("Medium", fMediumSizeXY / 2, fMediumSizeXY / 2, fMediumSizeZ / 2);
0118
0119   fLogicMedium = new G4LogicalVolume(fSolidMedium, fMediumMaterial, "Medium");
0120
0121   fPhysiMedium = new G4PVPlacement(nullptr,
0122                                    G4ThreeVector(0, 0, 0),
0123                                    "Medium",
0124                                    fLogicMedium,
0125                                    fPhysiWorld,
0126                                    false,
0127                                    0);
0128
0129   // ************
0130   // Cell phantom
0131   // ************
0132
0133   // The cell phantom is placed in the middle of the parent volume (fLogicMedium here)
0134
0135   fPhantomParam = new CellParameterisation
0136     (fPhantomFileName, fRedMaterial, fGreenMaterial, fBlueMaterial, fShiftX, fShiftY, fShiftZ);
0137
0138   fSolidPhantom = new G4Box("Phantom",
0139                             fPhantomParam->GetPixelSizeX() / 2,
0140                             fPhantomParam->GetPixelSizeY() / 2,
0141                             fPhantomParam->GetPixelSizeZ() / 2);
0142
0143   fLogicPhantom = new G4LogicalVolume(fSolidPhantom,
0144                                       fPhantomMaterial, // material is not relevant,
0145                                                         // it will be changed by the
0146                                                         // ComputeMaterial method
0147                                                         // of the CellParameterisation
0148                                       "Phantom",
0149                                       nullptr,
0150                                       nullptr,
0151                                       nullptr);
0152
0153   fPhysiPhantom = new G4PVParameterised(
0154       "Phantom",         // name
0155       fLogicPhantom,     // logical volume
0156       fLogicMedium,      // mother logical volume
0157       kUndefined,        // kUndefined: three-dimensional optimization
0158       fPhantomParam->GetPhantomTotalPixels(), // number of voxels
0159       fPhantomParam,     // the parametrisation
0160       false);
0161
0162   G4cout << " #########################################################################" << G4endl;
0163   G4cout << "                               Phantom information                        " << G4endl;
0164   G4cout << " #########################################################################" << G4endl;
0165   G4cout << G4endl;
0166
0167   G4cout << " ==========> The phantom contains " << fPhantomParam->GetPhantomTotalPixels()
0168          << " voxels " << G4endl;
0169   G4cout << " ==========> Voxel size X (um) = " << fPhantomParam->GetPixelSizeX()/um << G4endl;
0170   G4cout << " ==========> Voxel size Y (um) = " << fPhantomParam->GetPixelSizeY()/um << G4endl;
0171   G4cout << " ==========> Voxel size Z (um) = " << fPhantomParam->GetPixelSizeZ()/um << G4endl;
0172   G4cout << G4endl;
0173
0174   G4cout << " ==========> Number of red voxels = "
0175          << fPhantomParam->GetRedTotalPixels() << G4endl;
0176   G4cout << " ==========> Number of green voxels = "
0177          << fPhantomParam->GetGreenTotalPixels() << G4endl;
0178   G4cout << " ==========> Number of blue voxels = "
0179          << fPhantomParam->GetBlueTotalPixels() << G4endl;
0180   G4cout << G4endl;
0181
0182   G4cout << " ==========> Tolal mass of red voxels (kg) = "
0183          << fPhantomParam->GetRedMass() / kg << G4endl;
0184   G4cout << " ==========> Tolal mass of green voxels (kg) = "
0185          << fPhantomParam->GetGreenMass() / kg << G4endl;
0186   G4cout << " ==========> Tolal mass of blue voxels (kg) = "
0187         << fPhantomParam->GetBlueMass() / kg << G4endl;
0188   G4cout << G4endl;
0189   G4cout << " #########################################################################" << G4endl;
0190   G4cout << G4endl;
0191
0192   // USER LIMITS ON STEP LENGTH
0193
0194   // fLogicWorld->SetUserLimits(new G4UserLimits(100 * mm));
0195   // fLogicPhantom->SetUserLimits(new G4UserLimits(0.5 * micrometer));
0196   // fLogicMedium->SetUserLimits(new G4UserLimits(1 * micrometer));
0197
0198   // Create a phantom G4Region and add logical volume
0199
0200   fPhantomRegion = new G4Region("phantomRegion");
0201
0202   G4ProductionCuts* cuts = new G4ProductionCuts();
0203
0204   G4double defCut = 1*nanometer;
0205   cuts->SetProductionCut(defCut,"gamma");
0206   cuts->SetProductionCut(defCut,"e-");
0207   cuts->SetProductionCut(defCut,"e+");
0208   cuts->SetProductionCut(defCut,"proton");
0209
0210   fPhantomRegion->SetProductionCuts(cuts);
0211   fPhantomRegion->AddRootLogicalVolume(fLogicMedium);
0212
0213   return fPhysiWorld;
0214 }
0215
0216 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0217
0218 void DetectorConstruction::SetTargetMaterial(const G4String& mat)
0219 {
0220   if (G4Material* material = G4NistManager::Instance()->FindOrBuildMaterial(mat))
0221   {
0222     if (material && mat != "G4_WATER")
0223     {
0224       fMediumMaterial = material;
0225       G4cout << " #########################################################################"
0226              << G4endl;
0227       G4cout << "                               Cell culture medium material               "
0228              << G4endl;
0229       G4cout << fMediumMaterial << G4endl;
0230       G4cout << " #########################################################################"
0231              << G4endl;
0232       G4cout << G4endl;
0233     }
0234   }
0235   else
0236   {
0237     G4cout << G4endl;
0238     G4cout << "WARNING: material \"" << mat << "\" doesn't exist in NIST elements/materials"
0239            << G4endl;
0240     G4cout << " table [located in $G4INSTALL/source/materials/src/G4NistMaterialBuilder.cc]"
0241            << G4endl;
0242     G4cout << G4endl;
0243   }
0244 }
0245
0246 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0247
0248 void DetectorConstruction::SetRedDensity(const G4double& value)
0249 {
0250   fDensityRed = value;
0251   if (fDensityRed != 1.0)
0252   {
0253     G4NistManager *man = G4NistManager::Instance();
0254     G4Material * H2O_red = man->BuildMaterialWithNewDensity("G4_WATER_red","G4_WATER",
0255                                                             fDensityRed);
0256     fRedMaterial = H2O_red;
0257   }
0258   else
0259   {
0260     G4NistManager *man = G4NistManager::Instance();
0261     fRedMaterial = man->FindOrBuildMaterial("G4_WATER");
0262   }
0263 }
0264
0265 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0266
0267 void DetectorConstruction::SetGreenDensity(const G4double& value)
0268 {
0269   fDensityGreen = value;
0270   if (fDensityGreen != 1.0)
0271   {
0272     G4NistManager *man = G4NistManager::Instance();
0273     G4Material * H2O_green = man->BuildMaterialWithNewDensity("G4_WATER_green","G4_WATER",
0274                                                               fDensityGreen);
0275     fGreenMaterial = H2O_green;
0276   }
0277   else
0278   {
0279     G4NistManager *man = G4NistManager::Instance();
0280     fGreenMaterial = man->FindOrBuildMaterial("G4_WATER");
0281   }
0282 }
0283
0284 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0285
0286 void DetectorConstruction::SetBlueDensity(const G4double& value)
0287 {
0288   fDensityBlue = value;
0289   if (fDensityBlue != 1.0)
0290   {
0291     G4NistManager *man = G4NistManager::Instance();
0292     G4Material * H2O_blue = man->BuildMaterialWithNewDensity("G4_WATER_blue","G4_WATER",
0293                                                              fDensityBlue);
0294     fBlueMaterial = H2O_blue;
0295   }
0296   else
0297   {
0298     G4NistManager *man = G4NistManager::Instance();
0299     fBlueMaterial = man->FindOrBuildMaterial("G4_WATER");
0300   }
0301 }
0302
0303 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0304
0305 void DetectorConstruction::SetShiftX(const G4double& value)
0306 {
0307   fShiftX = value;
0308   G4cout << "... setting phantom shift: X = " << fShiftX/um << " um" << G4endl;
0309 }
0310
0311 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0312
0313 void DetectorConstruction::SetShiftY(const G4double& value)
0314 {
0315   fShiftY = value;
0316   G4cout << "... setting phantom shift: Y = " << fShiftY/um << " um" << G4endl;
0317 }
0318
0319 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0320
0321 void DetectorConstruction::SetShiftZ(const G4double& value)
0322 {
0323   fShiftZ = value;
0324   G4cout << "... setting phantom shift: Y = " << fShiftZ/um << " um" << G4endl;
0325 }
0326
0327 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0328
0329 void DetectorConstruction::SetMediumSizeXY(const G4double& value)
0330 {
0331   fMediumSizeXY = value;
0332 }
0333
0334 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0335
0336 void DetectorConstruction::SetMediumSizeZ(const G4double& value)
0337 {
0338   fMediumSizeZ = value;
0339 }
0340
0341 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0342
0343 void DetectorConstruction::SetWorldSizeXY(const G4double& value)
0344 {
0345   fWorldSizeXY = value;
0346 }
0347
0348 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0349
0350 void DetectorConstruction::SetWorldSizeZ(const G4double& value)
0351 {
0352   fWorldSizeZ = value;
0353 }
0354
0355 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0356
0357 void DetectorConstruction::SetPhantomFileName(const G4String& phantomName)
0358 {
0359   fPhantomFileName = phantomName;
0360   G4cout << " #########################################################################"
0361          << G4endl;
0362   G4cout << "                          Loading cell phantom from file:  "
0363          << fPhantomFileName << G4endl;
0364   G4cout << " #########################################################################"
0365          << G4endl;
0366   G4cout << G4endl;
0367 }