| 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 /// \file medical/DICOM/src/DicomDetectorConstruction.cc
28 /// \brief Implementation of the DicomDetectorConstruction class
29 //
30
31 #include "DicomDetectorConstruction.hh"
32
33 #include "CLHEP/Units/SystemOfUnits.h"
34 #include "DicomHandler.hh"
35 #include "DicomPhantomZSliceHeader.hh"
36
37 #include "G4Box.hh"
38 #include "G4Element.hh"
39 #include "G4LogicalVolume.hh"
40 #include "G4Material.hh"
41 #include "G4NistManager.hh"
42 #include "G4PVPlacement.hh"
43 #include "G4PhysicalConstants.hh"
44 #include "G4UIcommand.hh"
45 #include "G4VPhysicalVolume.hh"
46 #include "globals.hh"
47
48 #ifdef G4_DCMTK
49 # include "DicomFileMgr.hh"
50 #endif
51 #include "G4VisAttributes.hh"
52
53 using CLHEP::cm3;
54 using CLHEP::g;
55 using CLHEP::m;
56 using CLHEP::mg;
57 using CLHEP::mole;
58 using CLHEP::perCent;
59
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
61 DicomDetectorConstruction::DicomDetectorConstruction()
62 : G4VUserDetectorConstruction(),
63 fAir(0),
64
65 fWorld_solid(0),
66 fWorld_logic(0),
67 fWorld_phys(0),
68
69 fContainer_solid(0),
70 fContainer_logic(0),
71 fContainer_phys(0),
72
73 fNoFiles(0),
74 fMateIDs(0),
75
76 fZSliceHeaderMerged(0),
77
78 fNoVoxelsX(0),
79 fNoVoxelsY(0),
80 fNoVoxelsZ(0),
81 fVoxelHalfDimX(0),
82 fVoxelHalfDimY(0),
83 fVoxelHalfDimZ(0),
84
85 fConstructed(false)
86 {}
87
88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
89 DicomDetectorConstruction::~DicomDetectorConstruction() {}
90
91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
92 G4VPhysicalVolume* DicomDetectorConstruction::Construct()
93 {
94 if (!fConstructed || fWorld_phys == 0) {
95 fConstructed = true;
96 InitialisationOfMaterials();
97
98 //----- Build world
99 G4double worldXDimension = 1. * m;
100 G4double worldYDimension = 1. * m;
101 G4double worldZDimension = 1. * m;
102
103 fWorld_solid = new G4Box("WorldSolid", worldXDimension, worldYDimension, worldZDimension);
104
105 fWorld_logic = new G4LogicalVolume(fWorld_solid, fAir, "WorldLogical", 0, 0, 0);
106
107 fWorld_phys = new G4PVPlacement(0, G4ThreeVector(0, 0, 0), "World", fWorld_logic, 0, false, 0);
108
109 fWorld_logic->SetVisAttributes(G4VisAttributes::GetInvisible());
110
111 #ifdef G4_DCMTK
112 ReadPhantomDataNew();
113 ConstructPhantomContainerNew();
114 #else
115 ReadPhantomData();
116 ConstructPhantomContainer();
117 #endif
118
119 ConstructPhantom();
120 }
121 return fWorld_phys;
122 }
123
124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........
125 void DicomDetectorConstruction::InitialisationOfMaterials()
126 {
127 // Creating elements :
128 G4double z, a, density;
129 G4String name, symbol;
130
131 G4Element* elC = new G4Element(name = "Carbon", symbol = "C", z = 6.0, a = 12.011 * g / mole);
132 G4Element* elH = new G4Element(name = "Hydrogen", symbol = "H", z = 1.0, a = 1.008 * g / mole);
133 G4Element* elN = new G4Element(name = "Nitrogen", symbol = "N", z = 7.0, a = 14.007 * g / mole);
134 G4Element* elO = new G4Element(name = "Oxygen", symbol = "O", z = 8.0, a = 16.00 * g / mole);
135 G4Element* elNa =
136 new G4Element(name = "Sodium", symbol = "Na", z = 11.0, a = 22.98977 * g / mole);
137 G4Element* elMg =
138 new G4Element(name = "Magnesium", symbol = "Mg", z = 12.0, a = 24.3050 * g / mole);
139 G4Element* elP =
140 new G4Element(name = "Phosphorus", symbol = "P", z = 15.0, a = 30.973976 * g / mole);
141 G4Element* elS = new G4Element(name = "Sulfur", symbol = "S", z = 16.0, a = 32.065 * g / mole);
142 G4Element* elCl =
143 new G4Element(name = "Chlorine", symbol = "Cl", z = 17.0, a = 35.453 * g / mole);
144 G4Element* elK =
145 new G4Element(name = "Potassium", symbol = "K", z = 19.0, a = 30.0983 * g / mole);
146
147 G4Element* elFe = new G4Element(name = "Iron", symbol = "Fe", z = 26, a = 56.845 * g / mole);
148
149 G4Element* elCa = new G4Element(name = "Calcium", symbol = "Ca", z = 20.0, a = 40.078 * g / mole);
150
151 G4Element* elZn = new G4Element(name = "Zinc", symbol = "Zn", z = 30.0, a = 65.382 * g / mole);
152
153 // Creating Materials :
154 G4int numberofElements;
155
156 // Air
157 fAir = new G4Material("Air", 1.290 * mg / cm3, numberofElements = 2);
158 fAir->AddElement(elN, 0.7);
159 fAir->AddElement(elO, 0.3);
160
161 // Soft tissue (ICRP - NIST)
162 G4Material* softTissue = new G4Material("SoftTissue", 1.00 * g / cm3, numberofElements = 13);
163 softTissue->AddElement(elH, 10.4472 * perCent);
164 softTissue->AddElement(elC, 23.219 * perCent);
165 softTissue->AddElement(elN, 2.488 * perCent);
166 softTissue->AddElement(elO, 63.0238 * perCent);
167 softTissue->AddElement(elNa, 0.113 * perCent);
168 softTissue->AddElement(elMg, 0.0113 * perCent);
169 softTissue->AddElement(elP, 0.113 * perCent);
170 softTissue->AddElement(elS, 0.199 * perCent);
171 softTissue->AddElement(elCl, 0.134 * perCent);
172 softTissue->AddElement(elK, 0.199 * perCent);
173 softTissue->AddElement(elCa, 0.023 * perCent);
174 softTissue->AddElement(elFe, 0.005 * perCent);
175 softTissue->AddElement(elZn, 0.003 * perCent);
176
177 // Lung Inhale
178 G4Material* lunginhale =
179 new G4Material("LungInhale", density = 0.217 * g / cm3, numberofElements = 9);
180 lunginhale->AddElement(elH, 0.103);
181 lunginhale->AddElement(elC, 0.105);
182 lunginhale->AddElement(elN, 0.031);
183 lunginhale->AddElement(elO, 0.749);
184 lunginhale->AddElement(elNa, 0.002);
185 lunginhale->AddElement(elP, 0.002);
186 lunginhale->AddElement(elS, 0.003);
187 lunginhale->AddElement(elCl, 0.002);
188 lunginhale->AddElement(elK, 0.003);
189
190 // Lung exhale
191 G4Material* lungexhale =
192 new G4Material("LungExhale", density = 0.508 * g / cm3, numberofElements = 9);
193 lungexhale->AddElement(elH, 0.103);
194 lungexhale->AddElement(elC, 0.105);
195 lungexhale->AddElement(elN, 0.031);
196 lungexhale->AddElement(elO, 0.749);
197 lungexhale->AddElement(elNa, 0.002);
198 lungexhale->AddElement(elP, 0.002);
199 lungexhale->AddElement(elS, 0.003);
200 lungexhale->AddElement(elCl, 0.002);
201 lungexhale->AddElement(elK, 0.003);
202
203 // Adipose tissue
204 G4Material* adiposeTissue =
205 new G4Material("AdiposeTissue", density = 0.967 * g / cm3, numberofElements = 7);
206 adiposeTissue->AddElement(elH, 0.114);
207 adiposeTissue->AddElement(elC, 0.598);
208 adiposeTissue->AddElement(elN, 0.007);
209 adiposeTissue->AddElement(elO, 0.278);
210 adiposeTissue->AddElement(elNa, 0.001);
211 adiposeTissue->AddElement(elS, 0.001);
212 adiposeTissue->AddElement(elCl, 0.001);
213
214 // Brain (ICRP - NIST)
215 G4Material* brainTissue = new G4Material("BrainTissue", 1.03 * g / cm3, numberofElements = 13);
216 brainTissue->AddElement(elH, 11.0667 * perCent);
217 brainTissue->AddElement(elC, 12.542 * perCent);
218 brainTissue->AddElement(elN, 1.328 * perCent);
219 brainTissue->AddElement(elO, 73.7723 * perCent);
220 brainTissue->AddElement(elNa, 0.1840 * perCent);
221 brainTissue->AddElement(elMg, 0.015 * perCent);
222 brainTissue->AddElement(elP, 0.356 * perCent);
223 brainTissue->AddElement(elS, 0.177 * perCent);
224 brainTissue->AddElement(elCl, 0.236 * perCent);
225 brainTissue->AddElement(elK, 0.31 * perCent);
226 brainTissue->AddElement(elCa, 0.009 * perCent);
227 brainTissue->AddElement(elFe, 0.005 * perCent);
228 brainTissue->AddElement(elZn, 0.001 * perCent);
229
230 // Breast
231 G4Material* breast = new G4Material("Breast", density = 0.990 * g / cm3, numberofElements = 8);
232 breast->AddElement(elH, 0.109);
233 breast->AddElement(elC, 0.506);
234 breast->AddElement(elN, 0.023);
235 breast->AddElement(elO, 0.358);
236 breast->AddElement(elNa, 0.001);
237 breast->AddElement(elP, 0.001);
238 breast->AddElement(elS, 0.001);
239 breast->AddElement(elCl, 0.001);
240
241 // Spinal Disc
242 G4Material* spinalDisc = new G4Material("SpinalDisc", 1.10 * g / cm3, numberofElements = 8);
243 spinalDisc->AddElement(elH, 9.60 * perCent);
244 spinalDisc->AddElement(elC, 9.90 * perCent);
245 spinalDisc->AddElement(elN, 2.20 * perCent);
246 spinalDisc->AddElement(elO, 74.40 * perCent);
247 spinalDisc->AddElement(elNa, 0.50 * perCent);
248 spinalDisc->AddElement(elP, 2.20 * perCent);
249 spinalDisc->AddElement(elS, 0.90 * perCent);
250 spinalDisc->AddElement(elCl, 0.30 * perCent);
251
252 // Water
253 G4Material* water = new G4Material("Water", density = 1.0 * g / cm3, numberofElements = 2);
254 water->AddElement(elH, 0.112);
255 water->AddElement(elO, 0.888);
256
257 // Muscle
258 G4Material* muscle = new G4Material("Muscle", density = 1.061 * g / cm3, numberofElements = 9);
259 muscle->AddElement(elH, 0.102);
260 muscle->AddElement(elC, 0.143);
261 muscle->AddElement(elN, 0.034);
262 muscle->AddElement(elO, 0.710);
263 muscle->AddElement(elNa, 0.001);
264 muscle->AddElement(elP, 0.002);
265 muscle->AddElement(elS, 0.003);
266 muscle->AddElement(elCl, 0.001);
267 muscle->AddElement(elK, 0.004);
268
269 // Liver
270 G4Material* liver = new G4Material("Liver", density = 1.071 * g / cm3, numberofElements = 9);
271 liver->AddElement(elH, 0.102);
272 liver->AddElement(elC, 0.139);
273 liver->AddElement(elN, 0.030);
274 liver->AddElement(elO, 0.716);
275 liver->AddElement(elNa, 0.002);
276 liver->AddElement(elP, 0.003);
277 liver->AddElement(elS, 0.003);
278 liver->AddElement(elCl, 0.002);
279 liver->AddElement(elK, 0.003);
280
281 // Tooth Dentin
282 G4Material* toothDentin = new G4Material("ToothDentin", 2.14 * g / cm3, numberofElements = 10);
283 toothDentin->AddElement(elH, 2.67 * perCent);
284 toothDentin->AddElement(elC, 12.77 * perCent);
285 toothDentin->AddElement(elN, 4.27 * perCent);
286 toothDentin->AddElement(elO, 40.40 * perCent);
287 toothDentin->AddElement(elNa, 0.65 * perCent);
288 toothDentin->AddElement(elMg, 0.59 * perCent);
289 toothDentin->AddElement(elP, 11.86 * perCent);
290 toothDentin->AddElement(elCl, 0.04 * perCent);
291 toothDentin->AddElement(elCa, 26.74 * perCent);
292 toothDentin->AddElement(elZn, 0.01 * perCent);
293
294 // Trabecular Bone
295 G4Material* trabecularBone =
296 new G4Material("TrabecularBone", density = 1.159 * g / cm3, numberofElements = 12);
297 trabecularBone->AddElement(elH, 0.085);
298 trabecularBone->AddElement(elC, 0.404);
299 trabecularBone->AddElement(elN, 0.058);
300 trabecularBone->AddElement(elO, 0.367);
301 trabecularBone->AddElement(elNa, 0.001);
302 trabecularBone->AddElement(elMg, 0.001);
303 trabecularBone->AddElement(elP, 0.034);
304 trabecularBone->AddElement(elS, 0.002);
305 trabecularBone->AddElement(elCl, 0.002);
306 trabecularBone->AddElement(elK, 0.001);
307 trabecularBone->AddElement(elCa, 0.044);
308 trabecularBone->AddElement(elFe, 0.001);
309
310 // Trabecular bone used in the DICOM Head
311
312 G4Material* trabecularBone_head =
313 new G4Material("TrabecularBone_HEAD", 1.18 * g / cm3, numberofElements = 12);
314 trabecularBone_head->AddElement(elH, 8.50 * perCent);
315 trabecularBone_head->AddElement(elC, 40.40 * perCent);
316 trabecularBone_head->AddElement(elN, 2.80 * perCent);
317 trabecularBone_head->AddElement(elO, 36.70 * perCent);
318 trabecularBone_head->AddElement(elNa, 0.10 * perCent);
319 trabecularBone_head->AddElement(elMg, 0.10 * perCent);
320 trabecularBone_head->AddElement(elP, 3.40 * perCent);
321 trabecularBone_head->AddElement(elS, 0.20 * perCent);
322 trabecularBone_head->AddElement(elCl, 0.20 * perCent);
323 trabecularBone_head->AddElement(elK, 0.10 * perCent);
324 trabecularBone_head->AddElement(elCa, 7.40 * perCent);
325 trabecularBone_head->AddElement(elFe, 0.10 * perCent);
326
327 // Dense Bone
328 G4Material* denseBone =
329 new G4Material("DenseBone", density = 1.575 * g / cm3, numberofElements = 11);
330 denseBone->AddElement(elH, 0.056);
331 denseBone->AddElement(elC, 0.235);
332 denseBone->AddElement(elN, 0.050);
333 denseBone->AddElement(elO, 0.434);
334 denseBone->AddElement(elNa, 0.001);
335 denseBone->AddElement(elMg, 0.001);
336 denseBone->AddElement(elP, 0.072);
337 denseBone->AddElement(elS, 0.003);
338 denseBone->AddElement(elCl, 0.001);
339 denseBone->AddElement(elK, 0.001);
340 denseBone->AddElement(elCa, 0.146);
341
342 // Cortical Bone (ICRP - NIST)
343 G4Material* corticalBone = new G4Material("CorticalBone", 1.85 * g / cm3, numberofElements = 9);
344 corticalBone->AddElement(elH, 4.7234 * perCent);
345 corticalBone->AddElement(elC, 14.4330 * perCent);
346 corticalBone->AddElement(elN, 4.199 * perCent);
347 corticalBone->AddElement(elO, 44.6096 * perCent);
348 corticalBone->AddElement(elMg, 0.22 * perCent);
349 corticalBone->AddElement(elP, 10.497 * perCent);
350 corticalBone->AddElement(elS, 0.315 * perCent);
351 corticalBone->AddElement(elCa, 20.993 * perCent);
352 corticalBone->AddElement(elZn, 0.01 * perCent);
353
354 // Tooth enamel
355 G4Material* toothEnamel = new G4Material("ToothEnamel", 2.89 * g / cm3, numberofElements = 10);
356 toothEnamel->AddElement(elH, 0.95 * perCent);
357 toothEnamel->AddElement(elC, 1.11 * perCent);
358 toothEnamel->AddElement(elN, 0.23 * perCent);
359 toothEnamel->AddElement(elO, 41.66 * perCent);
360 toothEnamel->AddElement(elNa, 0.79 * perCent);
361 toothEnamel->AddElement(elMg, 0.23 * perCent);
362 toothEnamel->AddElement(elP, 18.71 * perCent);
363 toothEnamel->AddElement(elCl, 0.34 * perCent);
364 toothEnamel->AddElement(elCa, 35.97 * perCent);
365 toothEnamel->AddElement(elZn, 0.02 * perCent);
366
367 #ifdef DICOM_USE_HEAD
368 //----- Put the materials in a vector HEAD PHANTOM
369 fOriginalMaterials.push_back(fAir); // 0.00129 g/cm3
370 fOriginalMaterials.push_back(softTissue); // 1.055 g/cm3
371 fOriginalMaterials.push_back(brainTissue); // 1.07 g/cm3
372 fOriginalMaterials.push_back(spinalDisc); // 1.10 g/cm3
373 fOriginalMaterials.push_back(trabecularBone_head); // 1.13 g/cm3
374 fOriginalMaterials.push_back(toothDentin); // 1.66 g/cm3
375 fOriginalMaterials.push_back(corticalBone); // 1.75 g/cm3
376 fOriginalMaterials.push_back(toothEnamel); // 2.04 g/cm3
377 G4cout << "The materials of the DICOM Head have been used" << G4endl;
378 #else
379 fOriginalMaterials.push_back(fAir); // rho = 0.00129
380 fOriginalMaterials.push_back(lunginhale); // rho = 0.217
381 fOriginalMaterials.push_back(lungexhale); // rho = 0.508
382 fOriginalMaterials.push_back(adiposeTissue); // rho = 0.967
383 fOriginalMaterials.push_back(breast); // rho = 0.990
384 fOriginalMaterials.push_back(water); // rho = 1.018
385 fOriginalMaterials.push_back(muscle); // rho = 1.061
386 fOriginalMaterials.push_back(liver); // rho = 1.071
387 fOriginalMaterials.push_back(trabecularBone); // rho = 1.159 - HEAD PHANTOM
388 fOriginalMaterials.push_back(denseBone); // rho = 1.575
389 G4cout << "Default materials of the DICOM Extended examples have been used" << G4endl;
390 #endif
391 }
392
393 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
394 void DicomDetectorConstruction::ReadPhantomDataNew()
395 {
396 #ifdef G4_DCMTK
397 G4String fileName = DicomFileMgr::GetInstance()->GetFileOutName();
398
399 std::ifstream fin(fileName);
400 std::vector<G4String> wl;
401 G4int nMaterials;
402 fin >> nMaterials;
403 G4String mateName;
404 G4int nmate;
405 for (G4int ii = 0; ii < nMaterials; ++ii) {
406 fin >> nmate;
407 fin >> mateName;
408 if (mateName[0] == '"' && mateName[mateName.length() - 1] == '"') {
409 mateName = mateName.substr(1, mateName.length() - 2);
410 }
411 G4cout << "GmReadPhantomG4Geometry::ReadPhantomData reading nmate " << ii << " = " << nmate
412 << " mate " << mateName << G4endl;
413 if (ii != nmate)
414 G4Exception("GmReadPhantomG4Geometry::ReadPhantomData", "Wrong argument",
415 FatalErrorInArgument,
416 "Material number should be in increasing order:wrong material number");
417
418 G4Material* mate = 0;
419 const G4MaterialTable* matTab = G4Material::GetMaterialTable();
420 for (auto matite = matTab->cbegin(); matite != matTab->cend(); ++matite) {
421 if ((*matite)->GetName() == mateName) {
422 mate = *matite;
423 }
424 }
425 if (mate == 0) {
426 mate = G4NistManager::Instance()->FindOrBuildMaterial(mateName);
427 }
428 if (!mate)
429 G4Exception("GmReadPhantomG4Geometry::ReadPhantomData", "Wrong argument",
430 FatalErrorInArgument, ("Material not found" + mateName).c_str());
431 fPhantomMaterialsOriginal[nmate] = mate;
432 }
433
434 fin >> fNoVoxelsX >> fNoVoxelsY >> fNoVoxelsZ;
435 G4cout << "GmReadPhantomG4Geometry::ReadPhantomData fNoVoxels X/Y/Z " << fNoVoxelsX << " "
436 << fNoVoxelsY << " " << fNoVoxelsZ << G4endl;
437 fin >> fMinX >> fMaxX;
438 fin >> fMinY >> fMaxY;
439 fin >> fMinZ >> fMaxZ;
440 fVoxelHalfDimX = (fMaxX - fMinX) / fNoVoxelsX / 2.;
441 fVoxelHalfDimY = (fMaxY - fMinY) / fNoVoxelsY / 2.;
442 fVoxelHalfDimZ = (fMaxZ - fMinZ) / fNoVoxelsZ / 2.;
443 # ifdef G4VERBOSE
444 G4cout << " Extension in X " << fMinX << " " << fMaxX << G4endl << " Extension in Y " << fMinY
445 << " " << fMaxY << G4endl << " Extension in Z " << fMinZ << " " << fMaxZ << G4endl;
446 # endif
447
448 fMateIDs = new size_t[fNoVoxelsX * fNoVoxelsY * fNoVoxelsZ];
449 for (G4int iz = 0; iz < fNoVoxelsZ; ++iz) {
450 for (G4int iy = 0; iy < fNoVoxelsY; ++iy) {
451 for (G4int ix = 0; ix < fNoVoxelsX; ++ix) {
452 G4int mateID;
453 fin >> mateID;
454 G4int nnew = ix + (iy)*fNoVoxelsX + (iz)*fNoVoxelsX * fNoVoxelsY;
455 if (mateID < 0 || mateID >= nMaterials) {
456 G4Exception("GmReadPhantomG4Geometry::ReadPhantomData", "Wrong index in phantom file",
457 FatalException,
458 G4String("It should be between 0 and "
459 + G4UIcommand::ConvertToString(nMaterials - 1) + ", while it is "
460 + G4UIcommand::ConvertToString(mateID))
461 .c_str());
462 }
463 fMateIDs[nnew] = mateID;
464 }
465 }
466 }
467
468 ReadVoxelDensities(fin);
469
470 fin.close();
471 #endif
472 }
473
474 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
475 void DicomDetectorConstruction::ReadVoxelDensities(std::ifstream& fin)
476 {
477 G4String stemp;
478 std::map<G4int, std::pair<G4double, G4double>> densiMinMax;
479 std::map<G4int, std::pair<G4double, G4double>>::iterator mpite;
480 for (G4int ii = 0; ii < G4int(fPhantomMaterialsOriginal.size()); ++ii) {
481 densiMinMax[ii] = std::pair<G4double, G4double>(DBL_MAX, -DBL_MAX);
482 }
483
484 char* part = std::getenv("DICOM_CHANGE_MATERIAL_DENSITY");
485 G4double densityDiff = -1.;
486 if (part) densityDiff = G4UIcommand::ConvertToDouble(part);
487
488 std::map<G4int, G4double> densityDiffs;
489 for (G4int ii = 0; ii < G4int(fPhantomMaterialsOriginal.size()); ++ii) {
490 densityDiffs[ii] = densityDiff; // currently all materials with same step
491 }
492 // densityDiffs[0] = 0.0001; //air
493
494 //--- Calculate the average material density for each material/density bin
495 std::map<std::pair<G4Material*, G4int>, matInfo*> newMateDens;
496
497 //---- Read the material densities
498 G4double dens;
499 for (G4int iz = 0; iz < fNoVoxelsZ; ++iz) {
500 for (G4int iy = 0; iy < fNoVoxelsY; ++iy) {
501 for (G4int ix = 0; ix < fNoVoxelsX; ++ix) {
502 fin >> dens;
503 G4int copyNo = ix + (iy)*fNoVoxelsX + (iz)*fNoVoxelsX * fNoVoxelsY;
504
505 if (densityDiff != -1.) continue;
506
507 //--- store the minimum and maximum density for each material
508 mpite = densiMinMax.find(G4int(fMateIDs[copyNo]));
509 if (dens < (*mpite).second.first) (*mpite).second.first = dens;
510 if (dens > (*mpite).second.second) (*mpite).second.second = dens;
511 //--- Get material from original list of material in file
512 G4int mateID = G4int(fMateIDs[copyNo]);
513 std::map<G4int, G4Material*>::const_iterator imite = fPhantomMaterialsOriginal.find(mateID);
514
515 //--- Check if density is equal to the original material density
516 if (std::fabs(dens - (*imite).second->GetDensity() / CLHEP::g * CLHEP::cm3) < 1.e-9)
517 continue;
518
519 //--- Build material name with fPhantomMaterialsOriginal name+density
520 G4int densityBin = (G4int(dens / densityDiffs[mateID]));
521
522 G4String mateName = (*imite).second->GetName() + G4UIcommand::ConvertToString(densityBin);
523 //--- Look if it is the first voxel with this material/densityBin
524 std::pair<G4Material*, G4int> matdens((*imite).second, densityBin);
525
526 auto mppite = newMateDens.find(matdens);
527 if (mppite != newMateDens.cend()) {
528 matInfo* mi = (*mppite).second;
529 mi->fSumdens += dens;
530 mi->fNvoxels++;
531 fMateIDs[copyNo] = fPhantomMaterialsOriginal.size() - 1 + mi->fId;
532 }
533 else {
534 matInfo* mi = new matInfo;
535 mi->fSumdens = dens;
536 mi->fNvoxels = 1;
537 mi->fId = G4int(newMateDens.size() + 1);
538 newMateDens[matdens] = mi;
539 fMateIDs[copyNo] = fPhantomMaterialsOriginal.size() - 1 + mi->fId;
540 }
541 }
542 }
543 }
544
545 if (densityDiff != -1.) {
546 for (mpite = densiMinMax.begin(); mpite != densiMinMax.end(); ++mpite) {
547 #ifdef G4VERBOSE
548 G4cout << "DicomDetectorConstruction::ReadVoxelDensities"
549 << " ORIG MATERIALS DENSITY " << (*mpite).first << " MIN " << (*mpite).second.first
550 << " MAX " << (*mpite).second.second << G4endl;
551 #endif
552 }
553 }
554
555 //----- Build the list of phantom materials that go to Parameterisation
556 //--- Add original materials
557 for (auto mimite = fPhantomMaterialsOriginal.cbegin(); mimite != fPhantomMaterialsOriginal.cend();
558 ++mimite)
559 {
560 fMaterials.push_back((*mimite).second);
561 }
562 //
563 //---- Build and add new materials
564 for (auto mppite = newMateDens.cbegin(); mppite != newMateDens.cend(); ++mppite) {
565 G4double averdens = (*mppite).second->fSumdens / (*mppite).second->fNvoxels;
566 G4double saverdens = G4int(1000.001 * averdens) / 1000.;
567 #ifndef G4VERBOSE
568 G4cout << "DicomDetectorConstruction::ReadVoxelDensities AVER DENS " << averdens << " -> "
569 << saverdens << " -> " << G4int(1000 * averdens) << " " << G4int(1000 * averdens) / 1000
570 << " " << G4int(1000 * averdens) / 1000. << G4endl;
571 #endif
572
573 G4String mateName =
574 ((*mppite).first).first->GetName() + "_" + G4UIcommand::ConvertToString(saverdens);
575 fMaterials.push_back(
576 BuildMaterialWithChangingDensity((*mppite).first.first, G4float(averdens), mateName));
577 }
578 }
579
580 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.......
581 void DicomDetectorConstruction::ReadPhantomData()
582 {
583 G4String dataFile = DicomHandler::GetDicomDataFile();
584 std::ifstream finDF(dataFile.c_str());
585 G4String fname;
586
587 if (finDF.good() != 1) {
588 G4String descript = "Problem reading data file: " + dataFile;
589 G4Exception(" DicomDetectorConstruction::ReadPhantomData", " ", FatalException, descript);
590 }
591
592 G4int compression;
593 finDF >> compression; // not used here
594 finDF >> fNoFiles;
595
596 for (G4int i = 0; i < fNoFiles; ++i) {
597 finDF >> fname;
598
599 //--- Read one data file
600 fname += ".g4dcm";
601
602 ReadPhantomDataFile(fname);
603 }
604
605 //----- Merge data headers
606 MergeZSliceHeaders();
607 finDF.close();
608 }
609
610 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........
611 void DicomDetectorConstruction::ReadPhantomDataFile(const G4String& fname)
612 {
613 G4cout << " DicomDetectorConstruction::ReadPhantomDataFile opening file " << fname
614 << G4endl; // GDEB
615
616 #ifdef G4VERBOSE
617 G4cout << " DicomDetectorConstruction::ReadPhantomDataFile opening file " << fname << G4endl;
618 #endif
619
620 std::ifstream fin(fname.c_str(), std::ios_base::in);
621 if (!fin.is_open()) {
622 G4Exception("DicomDetectorConstruction::ReadPhantomDataFile", "", FatalErrorInArgument,
623 G4String("File not found " + fname).c_str());
624 }
625 //----- Define density differences (maximum density difference to create
626 // a new material)
627 char* part = std::getenv("DICOM_CHANGE_MATERIAL_DENSITY");
628 G4double densityDiff = -1.;
629 if (part) densityDiff = G4UIcommand::ConvertToDouble(part);
630 if (densityDiff != -1.) {
631 for (unsigned int ii = 0; ii < fOriginalMaterials.size(); ++ii) {
632 fDensityDiffs[ii] = densityDiff; // currently all materials with
633 // same difference
634 }
635 }
636 else {
637 if (fMaterials.size() == 0) { // do it only for first slice
638 for (unsigned int ii = 0; ii < fOriginalMaterials.size(); ++ii) {
639 fMaterials.push_back(fOriginalMaterials[ii]);
640 }
641 }
642 }
643
644 //----- Read data header
645 DicomPhantomZSliceHeader* sliceHeader = new DicomPhantomZSliceHeader(fin);
646 fZSliceHeaders.push_back(sliceHeader);
647
648 //----- Read material indices
649 G4int nVoxels = sliceHeader->GetNoVoxels();
650
651 //--- If first slice, initiliaze fMateIDs
652 if (fZSliceHeaders.size() == 1) {
653 // fMateIDs = new unsigned int[fNoFiles*nVoxels];
654 fMateIDs = new size_t[fNoFiles * nVoxels];
655 }
656
657 unsigned int mateID;
658 // number of voxels from previously read slices
659 G4int voxelCopyNo = G4int((fZSliceHeaders.size() - 1) * nVoxels);
660 for (G4int ii = 0; ii < nVoxels; ++ii, voxelCopyNo++) {
661 fin >> mateID;
662 fMateIDs[voxelCopyNo] = mateID;
663 }
664
665 //----- Read material densities and build new materials if two voxels have
666 // same material but its density is in a different density interval
667 // (size of density intervals defined by densityDiff)
668 G4double density;
669 // number of voxels from previously read slices
670 voxelCopyNo = G4int((fZSliceHeaders.size() - 1) * nVoxels);
671 for (G4int ii = 0; ii < nVoxels; ++ii, voxelCopyNo++) {
672 fin >> density;
673
674 //-- Get material from list of original materials
675 mateID = unsigned(fMateIDs[voxelCopyNo]);
676 G4Material* mateOrig = fOriginalMaterials[mateID];
677
678 //-- Get density bin: middle point of the bin in which the current
679 // density is included
680 G4String newMateName = mateOrig->GetName();
681 G4float densityBin = 0.;
682 if (densityDiff != -1.) {
683 densityBin = G4float(fDensityDiffs[mateID]) * (G4int(density / fDensityDiffs[mateID]) + 0.5);
684 //-- Build the new material name
685 newMateName += G4UIcommand::ConvertToString(densityBin);
686 }
687
688 //-- Look if a material with this name is already created
689 // (because a previous voxel was already in this density bin)
690 unsigned int im;
691 for (im = 0; im < fMaterials.size(); ++im) {
692 if (fMaterials[im]->GetName() == newMateName) {
693 break;
694 }
695 }
696 //-- If material is already created use index of this material
697 if (im != fMaterials.size()) {
698 fMateIDs[voxelCopyNo] = im;
699 //-- else, create the material
700 }
701 else {
702 if (densityDiff != -1.) {
703 fMaterials.push_back(BuildMaterialWithChangingDensity(mateOrig, densityBin, newMateName));
704 fMateIDs[voxelCopyNo] = fMaterials.size() - 1;
705 }
706 else {
707 G4cerr << " im " << im << " < " << fMaterials.size() << " name " << newMateName << G4endl;
708 G4Exception("DicomDetectorConstruction::ReadPhantomDataFile", "", FatalErrorInArgument,
709 "Wrong index in material"); // it should never reach here
710 }
711 }
712 }
713 }
714
715 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
716 void DicomDetectorConstruction::MergeZSliceHeaders()
717 {
718 //----- Images must have the same dimension ...
719 fZSliceHeaderMerged = new DicomPhantomZSliceHeader(*fZSliceHeaders[0]);
720 for (unsigned int ii = 1; ii < fZSliceHeaders.size(); ++ii) {
721 *fZSliceHeaderMerged += *fZSliceHeaders[ii];
722 }
723 }
724
725 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
726 G4Material* DicomDetectorConstruction::BuildMaterialWithChangingDensity(const G4Material* origMate,
727 G4float density,
728 G4String newMateName)
729 {
730 //----- Copy original material, but with new density
731 G4int nelem = G4int(origMate->GetNumberOfElements());
732 G4Material* mate =
733 new G4Material(newMateName, density * g / cm3, nelem, kStateUndefined, STP_Temperature);
734
735 for (G4int ii = 0; ii < nelem; ++ii) {
736 G4double frac = origMate->GetFractionVector()[ii];
737 G4Element* elem = const_cast<G4Element*>(origMate->GetElement(ii));
738 mate->AddElement(elem, frac);
739 }
740
741 return mate;
742 }
743
744 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
745 void DicomDetectorConstruction::ConstructPhantomContainer()
746 {
747 //---- Extract number of voxels and voxel dimensions
748 fNoVoxelsX = fZSliceHeaderMerged->GetNoVoxelsX();
749 fNoVoxelsY = fZSliceHeaderMerged->GetNoVoxelsY();
750 fNoVoxelsZ = fZSliceHeaderMerged->GetNoVoxelsZ();
751
752 fVoxelHalfDimX = fZSliceHeaderMerged->GetVoxelHalfX();
753 fVoxelHalfDimY = fZSliceHeaderMerged->GetVoxelHalfY();
754 fVoxelHalfDimZ = fZSliceHeaderMerged->GetVoxelHalfZ();
755 #ifdef G4VERBOSE
756 G4cout << " fNoVoxelsX " << fNoVoxelsX << " fVoxelHalfDimX " << fVoxelHalfDimX << G4endl;
757 G4cout << " fNoVoxelsY " << fNoVoxelsY << " fVoxelHalfDimY " << fVoxelHalfDimY << G4endl;
758 G4cout << " fNoVoxelsZ " << fNoVoxelsZ << " fVoxelHalfDimZ " << fVoxelHalfDimZ << G4endl;
759 G4cout << " totalPixels " << fNoVoxelsX * fNoVoxelsY * fNoVoxelsZ << G4endl;
760 #endif
761
762 //----- Define the volume that contains all the voxels
763 fContainer_solid = new G4Box("phantomContainer", fNoVoxelsX * fVoxelHalfDimX,
764 fNoVoxelsY * fVoxelHalfDimY, fNoVoxelsZ * fVoxelHalfDimZ);
765 fContainer_logic =
766 new G4LogicalVolume(fContainer_solid,
767 // the material is not important, it will be fully filled by the voxels
768 fMaterials[0], "phantomContainer", 0, 0, 0);
769 //--- Place it on the world
770 G4double fOffsetX = (fZSliceHeaderMerged->GetMaxX() + fZSliceHeaderMerged->GetMinX()) / 2.;
771 G4double fOffsetY = (fZSliceHeaderMerged->GetMaxY() + fZSliceHeaderMerged->GetMinY()) / 2.;
772 G4double fOffsetZ = (fZSliceHeaderMerged->GetMaxZ() + fZSliceHeaderMerged->GetMinZ()) / 2.;
773 G4ThreeVector posCentreVoxels(fOffsetX, fOffsetY, fOffsetZ);
774 #ifdef G4VERBOSE
775 G4cout << " placing voxel container volume at " << posCentreVoxels << G4endl;
776 #endif
777 fContainer_phys = new G4PVPlacement(0, // rotation
778 posCentreVoxels,
779 fContainer_logic, // The logic volume
780 "phantomContainer", // Name
781 fWorld_logic, // Mother
782 false, // No op. bool.
783 1); // Copy number
784 }
785
786 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
787 void DicomDetectorConstruction::ConstructPhantomContainerNew()
788 {
789 #ifdef G4_DCMTK
790 //---- Extract number of voxels and voxel dimensions
791 # ifdef G4VERBOSE
792 G4cout << " fNoVoxelsX " << fNoVoxelsX << " fVoxelHalfDimX " << fVoxelHalfDimX << G4endl;
793 G4cout << " fNoVoxelsY " << fNoVoxelsY << " fVoxelHalfDimY " << fVoxelHalfDimY << G4endl;
794 G4cout << " fNoVoxelsZ " << fNoVoxelsZ << " fVoxelHalfDimZ " << fVoxelHalfDimZ << G4endl;
795 G4cout << " totalPixels " << fNoVoxelsX * fNoVoxelsY * fNoVoxelsZ << G4endl;
796 # endif
797
798 //----- Define the volume that contains all the voxels
799 fContainer_solid = new G4Box("phantomContainer", fNoVoxelsX * fVoxelHalfDimX,
800 fNoVoxelsY * fVoxelHalfDimY, fNoVoxelsZ * fVoxelHalfDimZ);
801 fContainer_logic =
802 new G4LogicalVolume(fContainer_solid,
803 // the material is not important, it will be fully filled by the voxels
804 fMaterials[0], "phantomContainer", 0, 0, 0);
805
806 G4ThreeVector posCentreVoxels((fMinX + fMaxX) / 2., (fMinY + fMaxY) / 2., (fMinZ + fMaxZ) / 2.);
807 # ifdef G4VERBOSE
808 G4cout << " placing voxel container volume at " << posCentreVoxels << G4endl;
809 # endif
810 fContainer_phys = new G4PVPlacement(0, // rotation
811 posCentreVoxels,
812 fContainer_logic, // The logic volume
813 "phantomContainer", // Name
814 fWorld_logic, // Mother
815 false, // No op. bool.
816 1); // Copy number
817 #endif
818 }
819
820 #include "G4MultiFunctionalDetector.hh"
821 #include "G4PSDoseDeposit.hh"
822 #include "G4PSDoseDeposit3D.hh"
823 #include "G4SDManager.hh"
824
825 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.
826 void DicomDetectorConstruction::SetScorer(G4LogicalVolume* voxel_logic)
827 {
828 #ifdef G4VERBOSE
829 G4cout << "\t SET SCORER : " << voxel_logic->GetName() << G4endl;
830 #endif
831
832 fScorers.insert(voxel_logic);
833 }
834
835 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
836
837 void DicomDetectorConstruction::ConstructSDandField()
838 {
839 #ifdef G4VERBOSE
840 G4cout << "\t CONSTRUCT SD AND FIELD" << G4endl;
841 #endif
842
843 // G4SDManager* SDman = G4SDManager::GetSDMpointer();
844
845 // SDman->SetVerboseLevel(1);
846
847 //
848 // Sensitive Detector Name
849 G4String concreteSDname = "phantomSD";
850 std::vector<G4String> scorer_names;
851 scorer_names.push_back(concreteSDname);
852 //------------------------
853 // MultiFunctionalDetector
854 //------------------------
855 //
856 // Define MultiFunctionalDetector with name.
857 // declare MFDet as a MultiFunctionalDetector scorer
858 G4MultiFunctionalDetector* MFDet = new G4MultiFunctionalDetector(concreteSDname);
859 G4SDManager::GetSDMpointer()->AddNewDetector(MFDet);
860 char* nest = std::getenv("DICOM_NESTED_PARAM");
861 if (nest && G4String(nest) == "1") {
862 G4VPrimitiveScorer* dosedep = new G4PSDoseDeposit3D(
863 "DoseDeposit", fNoVoxelsZ, fNoVoxelsY, fNoVoxelsX, 0, 2, 1); // nested param replica indexing
864 // - the last 3 arguments correspond to the replica depth for Z, Y and X respectively
865 MFDet->RegisterPrimitive(dosedep);
866 }
867 else {
868 G4VPrimitiveScorer* dosedep = new G4PSDoseDeposit("DoseDeposit"); // regular geometry
869 MFDet->RegisterPrimitive(dosedep);
870 }
871
872 for (auto ite = fScorers.cbegin(); ite != fScorers.cend(); ++ite) {
873 SetSensitiveDetector(*ite, MFDet);
874 }
875 }
876