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