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