| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // --------------------------------------------------------------------------------
27 // MONTE CARLO SIMULATION OF REALISTIC GEOMETRY FROM MICROSCOPES IMAGES
28 //
29 // Authors and contributors:
30 // P. Barberet, S. Incerti, N. H. Tran, L. Morelli
31 //
32 // University of Bordeaux, CNRS, LP2i, UMR5797, Gradignan, France
33 //
34 // If you use this code, please cite the following publication:
35 // P. Barberet et al.,
36 // "Monte-Carlo dosimetry on a realistic cell monolayer
37 // geometry exposed to alpha particles."
38 // Ph. Barberet et al 2012 Phys. Med. Biol. 57 2189
39 // doi: 110.1088/0031-9155/57/8/2189
40 // --------------------------------------------------------------------------------
41
42
43 #include "DetectorConstruction.hh"
44 #include "DetectorMessenger.hh"
45
46 #include "G4PhysicalConstants.hh"
47 #include "G4NistManager.hh"
48 #include "G4ProductionCuts.hh"
49
50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
51
52 DetectorConstruction::DetectorConstruction()
53 :G4VUserDetectorConstruction()
54 {
55 fDetectorMessenger = new DetectorMessenger(this);
56 }
57
58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
59
60 G4VPhysicalVolume *DetectorConstruction::Construct()
61 {
62 DefineMaterials();
63 return ConstructLine();
64 }
65
66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
67
68 void DetectorConstruction::DefineMaterials()
69 {
70 G4String name, symbol;
71
72 // Water and air are defined from NIST material database
73 G4NistManager *man = G4NistManager::Instance();
74
75 G4Material *H2O = man->FindOrBuildMaterial("G4_WATER");
76 G4Material *Air = man->FindOrBuildMaterial("G4_AIR");
77
78 fDefaultMaterial = Air;
79 fPhantomMaterial = H2O; // material is not relevant
80 // it will be changed by the ComputeMaterial
81 // method of the CellParameterisation
82
83 // Default materials
84 if (fMediumMaterial == nullptr) {fMediumMaterial = H2O;}
85 if (fRedMaterial == nullptr) {fRedMaterial = H2O;}
86 if (fGreenMaterial == nullptr) {fGreenMaterial = H2O;}
87 if (fBlueMaterial == nullptr) {fBlueMaterial = H2O;}
88 }
89
90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
91
92 G4VPhysicalVolume *DetectorConstruction::ConstructLine() {
93
94 //*************
95 // World volume
96 //*************
97
98 fSolidWorld = new G4Box("World", //its name
99 fWorldSizeXY / 2, fWorldSizeXY / 2, fWorldSizeZ / 2); //its size
100
101 fLogicWorld = new G4LogicalVolume(fSolidWorld, //its solid
102 fDefaultMaterial, //its material
103 "World"); //its name
104
105 fPhysiWorld = new G4PVPlacement(nullptr, //no rotation
106 G4ThreeVector(), //at (0,0,0)
107 "World", //its name
108 fLogicWorld, //its logical volume
109 nullptr, //its mother volume
110 false, //no boolean operation
111 0); //copy number
112
113 //********************
114 // Cell culture medium
115 //********************
116
117 fSolidMedium = new G4Box("Medium", fMediumSizeXY / 2, fMediumSizeXY / 2, fMediumSizeZ / 2);
118
119 fLogicMedium = new G4LogicalVolume(fSolidMedium, fMediumMaterial, "Medium");
120
121 fPhysiMedium = new G4PVPlacement(nullptr,
122 G4ThreeVector(0, 0, 0),
123 "Medium",
124 fLogicMedium,
125 fPhysiWorld,
126 false,
127 0);
128
129 // ************
130 // Cell phantom
131 // ************
132
133 // The cell phantom is placed in the middle of the parent volume (fLogicMedium here)
134
135 fPhantomParam = new CellParameterisation
136 (fPhantomFileName, fRedMaterial, fGreenMaterial, fBlueMaterial, fShiftX, fShiftY, fShiftZ);
137
138 fSolidPhantom = new G4Box("Phantom",
139 fPhantomParam->GetPixelSizeX() / 2,
140 fPhantomParam->GetPixelSizeY() / 2,
141 fPhantomParam->GetPixelSizeZ() / 2);
142
143 fLogicPhantom = new G4LogicalVolume(fSolidPhantom,
144 fPhantomMaterial, // material is not relevant,
145 // it will be changed by the
146 // ComputeMaterial method
147 // of the CellParameterisation
148 "Phantom",
149 nullptr,
150 nullptr,
151 nullptr);
152
153 fPhysiPhantom = new G4PVParameterised(
154 "Phantom", // name
155 fLogicPhantom, // logical volume
156 fLogicMedium, // mother logical volume
157 kUndefined, // kUndefined: three-dimensional optimization
158 fPhantomParam->GetPhantomTotalPixels(), // number of voxels
159 fPhantomParam, // the parametrisation
160 false);
161
162 G4cout << " #########################################################################" << G4endl;
163 G4cout << " Phantom information " << G4endl;
164 G4cout << " #########################################################################" << G4endl;
165 G4cout << G4endl;
166
167 G4cout << " ==========> The phantom contains " << fPhantomParam->GetPhantomTotalPixels()
168 << " voxels " << G4endl;
169 G4cout << " ==========> Voxel size X (um) = " << fPhantomParam->GetPixelSizeX()/um << G4endl;
170 G4cout << " ==========> Voxel size Y (um) = " << fPhantomParam->GetPixelSizeY()/um << G4endl;
171 G4cout << " ==========> Voxel size Z (um) = " << fPhantomParam->GetPixelSizeZ()/um << G4endl;
172 G4cout << G4endl;
173
174 G4cout << " ==========> Number of red voxels = "
175 << fPhantomParam->GetRedTotalPixels() << G4endl;
176 G4cout << " ==========> Number of green voxels = "
177 << fPhantomParam->GetGreenTotalPixels() << G4endl;
178 G4cout << " ==========> Number of blue voxels = "
179 << fPhantomParam->GetBlueTotalPixels() << G4endl;
180 G4cout << G4endl;
181
182 G4cout << " ==========> Tolal mass of red voxels (kg) = "
183 << fPhantomParam->GetRedMass() / kg << G4endl;
184 G4cout << " ==========> Tolal mass of green voxels (kg) = "
185 << fPhantomParam->GetGreenMass() / kg << G4endl;
186 G4cout << " ==========> Tolal mass of blue voxels (kg) = "
187 << fPhantomParam->GetBlueMass() / kg << G4endl;
188 G4cout << G4endl;
189 G4cout << " #########################################################################" << G4endl;
190 G4cout << G4endl;
191
192 // USER LIMITS ON STEP LENGTH
193
194 // fLogicWorld->SetUserLimits(new G4UserLimits(100 * mm));
195 // fLogicPhantom->SetUserLimits(new G4UserLimits(0.5 * micrometer));
196 // fLogicMedium->SetUserLimits(new G4UserLimits(1 * micrometer));
197
198 // Create a phantom G4Region and add logical volume
199
200 fPhantomRegion = new G4Region("phantomRegion");
201
202 G4ProductionCuts* cuts = new G4ProductionCuts();
203
204 G4double defCut = 1*nanometer;
205 cuts->SetProductionCut(defCut,"gamma");
206 cuts->SetProductionCut(defCut,"e-");
207 cuts->SetProductionCut(defCut,"e+");
208 cuts->SetProductionCut(defCut,"proton");
209
210 fPhantomRegion->SetProductionCuts(cuts);
211 fPhantomRegion->AddRootLogicalVolume(fLogicMedium);
212
213 return fPhysiWorld;
214 }
215
216 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
217
218 void DetectorConstruction::SetTargetMaterial(const G4String& mat)
219 {
220 if (G4Material* material = G4NistManager::Instance()->FindOrBuildMaterial(mat))
221 {
222 if (material && mat != "G4_WATER")
223 {
224 fMediumMaterial = material;
225 G4cout << " #########################################################################"
226 << G4endl;
227 G4cout << " Cell culture medium material "
228 << G4endl;
229 G4cout << fMediumMaterial << G4endl;
230 G4cout << " #########################################################################"
231 << G4endl;
232 G4cout << G4endl;
233 }
234 }
235 else
236 {
237 G4cout << G4endl;
238 G4cout << "WARNING: material \"" << mat << "\" doesn't exist in NIST elements/materials"
239 << G4endl;
240 G4cout << " table [located in $G4INSTALL/source/materials/src/G4NistMaterialBuilder.cc]"
241 << G4endl;
242 G4cout << G4endl;
243 }
244 }
245
246 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
247
248 void DetectorConstruction::SetRedDensity(const G4double& value)
249 {
250 fDensityRed = value;
251 if (fDensityRed != 1.0)
252 {
253 G4NistManager *man = G4NistManager::Instance();
254 G4Material * H2O_red = man->BuildMaterialWithNewDensity("G4_WATER_red","G4_WATER",
255 fDensityRed);
256 fRedMaterial = H2O_red;
257 }
258 else
259 {
260 G4NistManager *man = G4NistManager::Instance();
261 fRedMaterial = man->FindOrBuildMaterial("G4_WATER");
262 }
263 }
264
265 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
266
267 void DetectorConstruction::SetGreenDensity(const G4double& value)
268 {
269 fDensityGreen = value;
270 if (fDensityGreen != 1.0)
271 {
272 G4NistManager *man = G4NistManager::Instance();
273 G4Material * H2O_green = man->BuildMaterialWithNewDensity("G4_WATER_green","G4_WATER",
274 fDensityGreen);
275 fGreenMaterial = H2O_green;
276 }
277 else
278 {
279 G4NistManager *man = G4NistManager::Instance();
280 fGreenMaterial = man->FindOrBuildMaterial("G4_WATER");
281 }
282 }
283
284 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
285
286 void DetectorConstruction::SetBlueDensity(const G4double& value)
287 {
288 fDensityBlue = value;
289 if (fDensityBlue != 1.0)
290 {
291 G4NistManager *man = G4NistManager::Instance();
292 G4Material * H2O_blue = man->BuildMaterialWithNewDensity("G4_WATER_blue","G4_WATER",
293 fDensityBlue);
294 fBlueMaterial = H2O_blue;
295 }
296 else
297 {
298 G4NistManager *man = G4NistManager::Instance();
299 fBlueMaterial = man->FindOrBuildMaterial("G4_WATER");
300 }
301 }
302
303 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
304
305 void DetectorConstruction::SetShiftX(const G4double& value)
306 {
307 fShiftX = value;
308 G4cout << "... setting phantom shift: X = " << fShiftX/um << " um" << G4endl;
309 }
310
311 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
312
313 void DetectorConstruction::SetShiftY(const G4double& value)
314 {
315 fShiftY = value;
316 G4cout << "... setting phantom shift: Y = " << fShiftY/um << " um" << G4endl;
317 }
318
319 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
320
321 void DetectorConstruction::SetShiftZ(const G4double& value)
322 {
323 fShiftZ = value;
324 G4cout << "... setting phantom shift: Y = " << fShiftZ/um << " um" << G4endl;
325 }
326
327 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
328
329 void DetectorConstruction::SetMediumSizeXY(const G4double& value)
330 {
331 fMediumSizeXY = value;
332 }
333
334 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
335
336 void DetectorConstruction::SetMediumSizeZ(const G4double& value)
337 {
338 fMediumSizeZ = value;
339 }
340
341 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
342
343 void DetectorConstruction::SetWorldSizeXY(const G4double& value)
344 {
345 fWorldSizeXY = value;
346 }
347
348 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
349
350 void DetectorConstruction::SetWorldSizeZ(const G4double& value)
351 {
352 fWorldSizeZ = value;
353 }
354
355 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
356
357 void DetectorConstruction::SetPhantomFileName(const G4String& phantomName)
358 {
359 fPhantomFileName = phantomName;
360 G4cout << " #########################################################################"
361 << G4endl;
362 G4cout << " Loading cell phantom from file: "
363 << fPhantomFileName << G4endl;
364 G4cout << " #########################################################################"
365 << G4endl;
366 G4cout << G4endl;
367 }
368