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