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