Loading...
Searching...
No Matches
Public Member Functions | Protected Member Functions | Protected Attributes | List of all members
DicomDetectorConstruction Class Referenceabstract

#include <Doxymodules_medical.h>

Inheritance diagram for DicomDetectorConstruction:
G4VUserDetectorConstruction DicomNestedParamDetectorConstruction DicomPartialDetectorConstruction DicomRegularDetectorConstruction

Public Member Functions

 DicomDetectorConstruction ()
 
 ~DicomDetectorConstruction ()
 
virtual G4VPhysicalVolumeConstruct ()
 
G4int GetTotalVoxels () const
 

Protected Member Functions

void InitialisationOfMaterials ()
 
void ReadPhantomData ()
 
void ReadPhantomDataNew ()
 
void ReadVoxelDensities (std::ifstream &fin)
 
void ReadPhantomDataFile (const G4String &fname)
 
void MergeZSliceHeaders ()
 
G4MaterialBuildMaterialWithChangingDensity (const G4Material *origMate, G4float density, G4String newMateName)
 
void ConstructPhantomContainer ()
 
void ConstructPhantomContainerNew ()
 
virtual void ConstructPhantom ()=0
 
void SetScorer (G4LogicalVolume *voxel_logic)
 
virtual void ConstructSDandField ()
 

Protected Attributes

G4MaterialfAir
 
G4BoxfWorld_solid
 
G4LogicalVolumefWorld_logic
 
G4VPhysicalVolumefWorld_phys
 
G4BoxfContainer_solid
 
G4LogicalVolumefContainer_logic
 
G4VPhysicalVolumefContainer_phys
 
G4int fNoFiles
 
std::vector< G4Material * > fOriginalMaterials
 
std::vector< G4Material * > fMaterials
 
size_t * fMateIDs
 
std::map< G4int, G4double > fDensityDiffs
 
std::vector< DicomPhantomZSliceHeader * > fZSliceHeaders
 
DicomPhantomZSliceHeaderfZSliceHeaderMerged
 
G4int fNoVoxelsX
 
G4int fNoVoxelsY
 
G4int fNoVoxelsZ
 
G4double fVoxelHalfDimX
 
G4double fVoxelHalfDimY
 
G4double fVoxelHalfDimZ
 
G4double fMinX
 
G4double fMinY
 
G4double fMinZ
 
G4double fMaxX
 
G4double fMaxY
 
G4double fMaxZ
 
std::map< G4int, G4Material * > fPhantomMaterialsOriginal
 
DicomPhantomZSliceMergedfMergedSlices
 
std::set< G4LogicalVolume * > fScorers
 
G4bool fConstructed
 

Detailed Description

Definition at line 19 of file Doxymodules_medical.h.

Constructor & Destructor Documentation

◆ DicomDetectorConstruction()

DicomDetectorConstruction::DicomDetectorConstruction ( )

Definition at line 61 of file DicomDetectorConstruction.cc.

63 fAir(0),
64
65 fWorld_solid(0),
66 fWorld_logic(0),
67 fWorld_phys(0),
68
72
73 fNoFiles(0),
74 fMateIDs(0),
75
77
78 fNoVoxelsX(0),
79 fNoVoxelsY(0),
80 fNoVoxelsZ(0),
84
85 fConstructed(false)
86{
87}
DicomPhantomZSliceHeader * fZSliceHeaderMerged

◆ ~DicomDetectorConstruction()

DicomDetectorConstruction::~DicomDetectorConstruction ( )

Definition at line 90 of file DicomDetectorConstruction.cc.

91{
92}

Member Function Documentation

◆ Construct()

G4VPhysicalVolume * DicomDetectorConstruction::Construct ( )
virtual

Reimplemented in DicomPartialDetectorConstruction.

Definition at line 95 of file DicomDetectorConstruction.cc.

96{
97 if(!fConstructed || fWorld_phys == 0) {
98 fConstructed = true;
100
101 //----- Build world
102 G4double worldXDimension = 1.*m;
103 G4double worldYDimension = 1.*m;
104 G4double worldZDimension = 1.*m;
105
106 fWorld_solid = new G4Box( "WorldSolid",
107 worldXDimension,
108 worldYDimension,
109 worldZDimension );
110
112 fAir,
113 "WorldLogical",
114 0, 0, 0 );
115
116 fWorld_phys = new G4PVPlacement( 0,
117 G4ThreeVector(0,0,0),
118 "World",
120 0,
121 false,
122 0 );
123
124 fWorld_logic->SetVisAttributes(G4VisAttributes::GetInvisible());
125
126#ifdef G4_DCMTK
129#else
132#endif
133
135 }
136 return fWorld_phys;
137}
virtual void ConstructPhantom()=0

◆ GetTotalVoxels()

G4int DicomDetectorConstruction::GetTotalVoxels ( ) const
inline

Definition at line 158 of file DicomDetectorConstruction.hh.

159{
161}

◆ InitialisationOfMaterials()

void DicomDetectorConstruction::InitialisationOfMaterials ( )
protected

Definition at line 140 of file DicomDetectorConstruction.cc.

141{
142 // Creating elements :
143 G4double z, a, density;
144 G4String name, symbol;
145
146 G4Element* elC = new G4Element( name = "Carbon",
147 symbol = "C",
148 z = 6.0, a = 12.011 * g/mole );
149 G4Element* elH = new G4Element( name = "Hydrogen",
150 symbol = "H",
151 z = 1.0, a = 1.008 * g/mole );
152 G4Element* elN = new G4Element( name = "Nitrogen",
153 symbol = "N",
154 z = 7.0, a = 14.007 * g/mole );
155 G4Element* elO = new G4Element( name = "Oxygen",
156 symbol = "O",
157 z = 8.0, a = 16.00 * g/mole );
158 G4Element* elNa = new G4Element( name = "Sodium",
159 symbol = "Na",
160 z= 11.0, a = 22.98977* g/mole );
161 G4Element* elMg = new G4Element( name = "Magnesium",
162 symbol = "Mg",
163 z = 12.0, a = 24.3050* g/mole );
164 G4Element* elP = new G4Element( name = "Phosphorus",
165 symbol = "P",
166 z = 15.0, a = 30.973976* g/mole );
167 G4Element* elS = new G4Element( name = "Sulfur",
168 symbol = "S",
169 z = 16.0,a = 32.065* g/mole );
170 G4Element* elCl = new G4Element( name = "Chlorine",
171 symbol = "P",
172 z = 17.0, a = 35.453* g/mole );
173 G4Element* elK = new G4Element( name = "Potassium",
174 symbol = "P",
175 z = 19.0, a = 30.0983* g/mole );
176
177 G4Element* elFe = new G4Element( name = "Iron",
178 symbol = "Fe",
179 z = 26, a = 56.845* g/mole );
180
181 G4Element* elCa = new G4Element( name="Calcium",
182 symbol = "Ca",
183 z = 20.0, a = 40.078* g/mole );
184
185 G4Element* elZn = new G4Element( name = "Zinc",
186 symbol = "Zn",
187 z = 30.0,a = 65.382* g/mole );
188
189 // Creating Materials :
190 G4int numberofElements;
191
192 // Air
193 fAir = new G4Material( "Air",
194 1.290*mg/cm3,
195 numberofElements = 2 );
196 fAir->AddElement(elN, 0.7);
197 fAir->AddElement(elO, 0.3);
198
199
200 // Soft tissue (ICRP - NIST)
201 G4Material* softTissue = new G4Material ("SoftTissue", 1.00*g/cm3,
202 numberofElements = 13);
203 softTissue->AddElement(elH, 10.4472*perCent);
204 softTissue->AddElement(elC, 23.219*perCent);
205 softTissue->AddElement(elN, 2.488*perCent);
206 softTissue->AddElement(elO, 63.0238*perCent);
207 softTissue->AddElement(elNa, 0.113*perCent);
208 softTissue->AddElement(elMg, 0.0113*perCent);
209 softTissue->AddElement(elP, 0.113*perCent);
210 softTissue->AddElement(elS, 0.199*perCent);
211 softTissue->AddElement(elCl, 0.134*perCent);
212 softTissue->AddElement(elK, 0.199*perCent);
213 softTissue->AddElement(elCa, 0.023*perCent);
214 softTissue->AddElement(elFe, 0.005*perCent);
215 softTissue->AddElement(elZn, 0.003*perCent);
216
217 // Lung Inhale
218 G4Material* lunginhale = new G4Material( "LungInhale",
219 density = 0.217*g/cm3,
220 numberofElements = 9);
221 lunginhale->AddElement(elH,0.103);
222 lunginhale->AddElement(elC,0.105);
223 lunginhale->AddElement(elN,0.031);
224 lunginhale->AddElement(elO,0.749);
225 lunginhale->AddElement(elNa,0.002);
226 lunginhale->AddElement(elP,0.002);
227 lunginhale->AddElement(elS,0.003);
228 lunginhale->AddElement(elCl,0.002);
229 lunginhale->AddElement(elK,0.003);
230
231 // Lung exhale
232 G4Material* lungexhale = new G4Material( "LungExhale",
233 density = 0.508*g/cm3,
234 numberofElements = 9 );
235 lungexhale->AddElement(elH,0.103);
236 lungexhale->AddElement(elC,0.105);
237 lungexhale->AddElement(elN,0.031);
238 lungexhale->AddElement(elO,0.749);
239 lungexhale->AddElement(elNa,0.002);
240 lungexhale->AddElement(elP,0.002);
241 lungexhale->AddElement(elS,0.003);
242 lungexhale->AddElement(elCl,0.002);
243 lungexhale->AddElement(elK,0.003);
244
245 // Adipose tissue
246 G4Material* adiposeTissue = new G4Material( "AdiposeTissue",
247 density = 0.967*g/cm3,
248 numberofElements = 7);
249 adiposeTissue->AddElement(elH,0.114);
250 adiposeTissue->AddElement(elC,0.598);
251 adiposeTissue->AddElement(elN,0.007);
252 adiposeTissue->AddElement(elO,0.278);
253 adiposeTissue->AddElement(elNa,0.001);
254 adiposeTissue->AddElement(elS,0.001);
255 adiposeTissue->AddElement(elCl,0.001);
256
257 // Brain (ICRP - NIST)
258 G4Material* brainTissue = new G4Material ("BrainTissue", 1.03 * g/cm3,
259 numberofElements = 13);
260 brainTissue->AddElement(elH, 11.0667*perCent);
261 brainTissue->AddElement(elC, 12.542*perCent);
262 brainTissue->AddElement(elN, 1.328*perCent);
263 brainTissue->AddElement(elO, 73.7723*perCent);
264 brainTissue->AddElement(elNa, 0.1840*perCent);
265 brainTissue->AddElement(elMg, 0.015*perCent);
266 brainTissue->AddElement(elP, 0.356*perCent);
267 brainTissue->AddElement(elS, 0.177*perCent);
268 brainTissue->AddElement(elCl, 0.236*perCent);
269 brainTissue->AddElement(elK, 0.31*perCent);
270 brainTissue->AddElement(elCa, 0.009*perCent);
271 brainTissue->AddElement(elFe, 0.005*perCent);
272 brainTissue->AddElement(elZn, 0.001*perCent);
273
274
275 // Breast
276 G4Material* breast = new G4Material( "Breast",
277 density = 0.990*g/cm3,
278 numberofElements = 8 );
279 breast->AddElement(elH,0.109);
280 breast->AddElement(elC,0.506);
281 breast->AddElement(elN,0.023);
282 breast->AddElement(elO,0.358);
283 breast->AddElement(elNa,0.001);
284 breast->AddElement(elP,0.001);
285 breast->AddElement(elS,0.001);
286 breast->AddElement(elCl,0.001);
287
288 // Spinal Disc
289 G4Material* spinalDisc = new G4Material ("SpinalDisc", 1.10 * g/cm3,
290 numberofElements = 8);
291 spinalDisc->AddElement(elH, 9.60*perCent);
292 spinalDisc->AddElement(elC, 9.90*perCent);
293 spinalDisc->AddElement(elN, 2.20*perCent);
294 spinalDisc->AddElement(elO, 74.40*perCent);
295 spinalDisc->AddElement(elNa, 0.50*perCent);
296 spinalDisc->AddElement(elP, 2.20*perCent);
297 spinalDisc->AddElement(elS, 0.90*perCent);
298 spinalDisc->AddElement(elCl, 0.30*perCent);
299
300
301 // Water
302 G4Material* water = new G4Material( "Water",
303 density = 1.0*g/cm3,
304 numberofElements = 2 );
305 water->AddElement(elH,0.112);
306 water->AddElement(elO,0.888);
307
308 // Muscle
309 G4Material* muscle = new G4Material( "Muscle",
310 density = 1.061*g/cm3,
311 numberofElements = 9 );
312 muscle->AddElement(elH,0.102);
313 muscle->AddElement(elC,0.143);
314 muscle->AddElement(elN,0.034);
315 muscle->AddElement(elO,0.710);
316 muscle->AddElement(elNa,0.001);
317 muscle->AddElement(elP,0.002);
318 muscle->AddElement(elS,0.003);
319 muscle->AddElement(elCl,0.001);
320 muscle->AddElement(elK,0.004);
321
322 // Liver
323 G4Material* liver = new G4Material( "Liver",
324 density = 1.071*g/cm3,
325 numberofElements = 9);
326 liver->AddElement(elH,0.102);
327 liver->AddElement(elC,0.139);
328 liver->AddElement(elN,0.030);
329 liver->AddElement(elO,0.716);
330 liver->AddElement(elNa,0.002);
331 liver->AddElement(elP,0.003);
332 liver->AddElement(elS,0.003);
333 liver->AddElement(elCl,0.002);
334 liver->AddElement(elK,0.003);
335
336 // Tooth Dentin
337 G4Material* toothDentin = new G4Material ("ToothDentin", 2.14 * g/cm3,
338 numberofElements = 10);
339 toothDentin->AddElement(elH, 2.67*perCent);
340 toothDentin->AddElement(elC, 12.77*perCent);
341 toothDentin->AddElement(elN, 4.27*perCent);
342 toothDentin->AddElement(elO, 40.40*perCent);
343 toothDentin->AddElement(elNa, 0.65*perCent);
344 toothDentin->AddElement(elMg, 0.59*perCent);
345 toothDentin->AddElement(elP, 11.86*perCent);
346 toothDentin->AddElement(elCl, 0.04*perCent);
347 toothDentin->AddElement(elCa, 26.74*perCent);
348 toothDentin->AddElement(elZn, 0.01*perCent);
349
350
351 // Trabecular Bone
352 G4Material* trabecularBone = new G4Material("TrabecularBone",
353 density = 1.159*g/cm3,
354 numberofElements = 12 );
355 trabecularBone->AddElement(elH,0.085);
356 trabecularBone->AddElement(elC,0.404);
357 trabecularBone->AddElement(elN,0.058);
358 trabecularBone->AddElement(elO,0.367);
359 trabecularBone->AddElement(elNa,0.001);
360 trabecularBone->AddElement(elMg,0.001);
361 trabecularBone->AddElement(elP,0.034);
362 trabecularBone->AddElement(elS,0.002);
363 trabecularBone->AddElement(elCl,0.002);
364 trabecularBone->AddElement(elK,0.001);
365 trabecularBone->AddElement(elCa,0.044);
366 trabecularBone->AddElement(elFe,0.001);
367
368 // Trabecular bone used in the DICOM Head
369
370 G4Material* trabecularBone_head = new G4Material ("TrabecularBone_HEAD",
371 1.18 * g/cm3,
372 numberofElements = 12);
373 trabecularBone_head->AddElement(elH, 8.50*perCent);
374 trabecularBone_head->AddElement(elC, 40.40*perCent);
375 trabecularBone_head->AddElement(elN, 2.80*perCent);
376 trabecularBone_head->AddElement(elO, 36.70*perCent);
377 trabecularBone_head->AddElement(elNa, 0.10*perCent);
378 trabecularBone_head->AddElement(elMg, 0.10*perCent);
379 trabecularBone_head->AddElement(elP, 3.40*perCent);
380 trabecularBone_head->AddElement(elS, 0.20*perCent);
381 trabecularBone_head->AddElement(elCl, 0.20*perCent);
382 trabecularBone_head->AddElement(elK, 0.10*perCent);
383 trabecularBone_head->AddElement(elCa, 7.40*perCent);
384 trabecularBone_head->AddElement(elFe, 0.10*perCent);
385
386 // Dense Bone
387 G4Material* denseBone = new G4Material( "DenseBone",
388 density = 1.575*g/cm3,
389 numberofElements = 11 );
390 denseBone->AddElement(elH,0.056);
391 denseBone->AddElement(elC,0.235);
392 denseBone->AddElement(elN,0.050);
393 denseBone->AddElement(elO,0.434);
394 denseBone->AddElement(elNa,0.001);
395 denseBone->AddElement(elMg,0.001);
396 denseBone->AddElement(elP,0.072);
397 denseBone->AddElement(elS,0.003);
398 denseBone->AddElement(elCl,0.001);
399 denseBone->AddElement(elK,0.001);
400 denseBone->AddElement(elCa,0.146);
401
402 // Cortical Bone (ICRP - NIST)
403 G4Material* corticalBone = new G4Material ("CorticalBone", 1.85 * g/cm3,
404 numberofElements = 9);
405 corticalBone->AddElement(elH, 4.7234*perCent);
406 corticalBone->AddElement(elC, 14.4330*perCent);
407 corticalBone->AddElement(elN, 4.199*perCent);
408 corticalBone->AddElement(elO, 44.6096*perCent);
409 corticalBone->AddElement(elMg, 0.22*perCent);
410 corticalBone->AddElement(elP, 10.497*perCent);
411 corticalBone->AddElement(elS, 0.315*perCent);
412 corticalBone->AddElement(elCa, 20.993*perCent);
413 corticalBone->AddElement(elZn, 0.01*perCent);
414
415
416 // Tooth enamel
417 G4Material* toothEnamel = new G4Material ("ToothEnamel", 2.89 * g/cm3,
418 numberofElements = 10);
419 toothEnamel->AddElement(elH, 0.95*perCent);
420 toothEnamel->AddElement(elC, 1.11*perCent);
421 toothEnamel->AddElement(elN, 0.23*perCent);
422 toothEnamel->AddElement(elO,41.66*perCent);
423 toothEnamel->AddElement(elNa, 0.79*perCent);
424 toothEnamel->AddElement(elMg, 0.23*perCent);
425 toothEnamel->AddElement(elP, 18.71*perCent);
426 toothEnamel->AddElement(elCl, 0.34*perCent);
427 toothEnamel->AddElement(elCa, 35.97*perCent);
428 toothEnamel->AddElement(elZn, 0.02*perCent);
429
430#ifdef DICOM_USE_HEAD
431 //----- Put the materials in a vector HEAD PHANTOM
432 fOriginalMaterials.push_back(fAir); //0.00129 g/cm3
433 fOriginalMaterials.push_back(softTissue); // 1.055 g/cm3
434 fOriginalMaterials.push_back(brainTissue); // 1.07 g/cm3
435 fOriginalMaterials.push_back(spinalDisc); // 1.10 g/cm3
436 fOriginalMaterials.push_back(trabecularBone_head); // 1.13 g/cm3
437 fOriginalMaterials.push_back(toothDentin); // 1.66 g/cm3
438 fOriginalMaterials.push_back(corticalBone); // 1.75 g/cm3
439 fOriginalMaterials.push_back(toothEnamel); // 2.04 g/cm3
440 G4cout << "The materials of the DICOM Head have been used" << G4endl;
441#else
442 fOriginalMaterials.push_back(fAir); // rho = 0.00129
443 fOriginalMaterials.push_back(lunginhale); // rho = 0.217
444 fOriginalMaterials.push_back(lungexhale); // rho = 0.508
445 fOriginalMaterials.push_back(adiposeTissue); // rho = 0.967
446 fOriginalMaterials.push_back(breast ); // rho = 0.990
447 fOriginalMaterials.push_back(water); // rho = 1.018
448 fOriginalMaterials.push_back(muscle); // rho = 1.061
449 fOriginalMaterials.push_back(liver); // rho = 1.071
450 fOriginalMaterials.push_back(trabecularBone); // rho = 1.159 - HEAD PHANTOM
451 fOriginalMaterials.push_back(denseBone); // rho = 1.575
452 G4cout << "Default materials of the DICOM Extended examples have been used"
453 << G4endl;
454#endif
455}
std::vector< ExP01TrackerHit * > a
std::vector< G4Material * > fOriginalMaterials

◆ ReadPhantomData()

void DicomDetectorConstruction::ReadPhantomData ( )
protected

Definition at line 661 of file DicomDetectorConstruction.cc.

662{
664 std::ifstream finDF(dataFile.c_str());
665 G4String fname;
666
667if(finDF.good() != 1 )
668 {
669 G4String descript = "Problem reading data file: "+dataFile;
670 G4Exception(" DicomDetectorConstruction::ReadPhantomData"," ",
671 FatalException,descript);
672 }
673
674 G4int compression;
675 finDF >> compression; // not used here
676 finDF >> fNoFiles;
677
678 for(G4int i = 0; i < fNoFiles; ++i )
679 {
680
681 finDF >> fname;
682
683 //--- Read one data file
684 fname += ".g4dcm";
685
686 ReadPhantomDataFile(fname);
687 }
688
689 //----- Merge data headers
691 finDF.close();
692}
void ReadPhantomDataFile(const G4String &fname)
static G4String GetDicomDataFile()

◆ ReadPhantomDataNew()

void DicomDetectorConstruction::ReadPhantomDataNew ( )
protected

Definition at line 458 of file DicomDetectorConstruction.cc.

459{
460#ifdef G4_DCMTK
462
463 std::ifstream fin(fileName);
464 std::vector<G4String> wl;
465 G4int nMaterials;
466 fin >> nMaterials;
467 G4String mateName;
468 G4int nmate;
469 for( G4int ii = 0; ii < nMaterials; ++ii )
470 {
471 fin >> nmate;
472 fin >> mateName;
473 if( mateName[0] == '"' && mateName[mateName.length()-1] == '"' ) {
474 mateName = mateName.substr(1,mateName.length()-2);
475 }
476 G4cout << "GmReadPhantomG4Geometry::ReadPhantomData reading nmate "
477 << ii << " = " << nmate
478 << " mate " << mateName << G4endl;
479 if( ii != nmate )
480 G4Exception("GmReadPhantomG4Geometry::ReadPhantomData",
481 "Wrong argument", FatalErrorInArgument,
482 "Material number should be in increasing order:wrong material number");
483
484 G4Material* mate = 0;
485 const G4MaterialTable* matTab = G4Material::GetMaterialTable();
486 for( auto matite = matTab->cbegin(); matite != matTab->cend(); ++matite )
487 {
488 if( (*matite)->GetName() == mateName ) {
489 mate = *matite;
490 }
491 }
492 if( mate == 0 ) {
493 mate = G4NistManager::Instance()->FindOrBuildMaterial(mateName);
494 }
495 if( !mate ) G4Exception("GmReadPhantomG4Geometry::ReadPhantomData",
496 "Wrong argument",
497 FatalErrorInArgument,
498 ("Material not found" + mateName).c_str());
499 fPhantomMaterialsOriginal[nmate] = mate;
500 }
501
502 fin >> fNoVoxelsX >> fNoVoxelsY >> fNoVoxelsZ;
503 G4cout << "GmReadPhantomG4Geometry::ReadPhantomData fNoVoxels X/Y/Z "
504 << fNoVoxelsX << " "
505 << fNoVoxelsY << " " << fNoVoxelsZ << G4endl;
506 fin >> fMinX >> fMaxX;
507 fin >> fMinY >> fMaxY;
508 fin >> fMinZ >> fMaxZ;
512#ifdef G4VERBOSE
513 G4cout << " Extension in X " << fMinX << " " << fMaxX << G4endl
514 << " Extension in Y " << fMinY << " " << fMaxY << G4endl
515 << " Extension in Z " << fMinZ << " " << fMaxZ << G4endl;
516#endif
517
519 for( G4int iz = 0; iz < fNoVoxelsZ; ++iz ) {
520 for( G4int iy = 0; iy < fNoVoxelsY; ++iy ) {
521 for( G4int ix = 0; ix < fNoVoxelsX; ++ix ) {
522 G4int mateID;
523 fin >> mateID;
524 G4int nnew = ix + (iy)*fNoVoxelsX + (iz)*fNoVoxelsX*fNoVoxelsY;
525 if( mateID < 0 || mateID >= nMaterials ) {
526 G4Exception("GmReadPhantomG4Geometry::ReadPhantomData",
527 "Wrong index in phantom file",
528 FatalException,
529 G4String("It should be between 0 and "
530 + G4UIcommand::ConvertToString(nMaterials-1)
531 + ", while it is "
532 + G4UIcommand::ConvertToString(mateID)).c_str());
533 }
534 fMateIDs[nnew] = mateID;
535 }
536 }
537 }
538
539 ReadVoxelDensities( fin );
540
541 fin.close();
542#endif
543
544}
std::map< G4int, G4Material * > fPhantomMaterialsOriginal
void ReadVoxelDensities(std::ifstream &fin)
static DicomFileMgr * GetInstance()
G4String GetFileOutName() const

◆ ReadVoxelDensities()

void DicomDetectorConstruction::ReadVoxelDensities ( std::ifstream &  fin)
protected

Definition at line 547 of file DicomDetectorConstruction.cc.

548{
549 G4String stemp;
550 std::map<G4int, std::pair<G4double,G4double> > densiMinMax;
551 std::map<G4int, std::pair<G4double,G4double> >::iterator mpite;
552 for( G4int ii = 0; ii < G4int(fPhantomMaterialsOriginal.size()); ++ii )
553 {
554 densiMinMax[ii] = std::pair<G4double,G4double>(DBL_MAX,-DBL_MAX);
555 }
556
557 char* part = std::getenv( "DICOM_CHANGE_MATERIAL_DENSITY" );
558 G4double densityDiff = -1.;
559 if( part ) densityDiff = G4UIcommand::ConvertToDouble(part);
560
561 std::map<G4int,G4double> densityDiffs;
562 for( G4int ii = 0; ii < G4int(fPhantomMaterialsOriginal.size()); ++ii )
563 {
564 densityDiffs[ii] = densityDiff; //currently all materials with same step
565 }
566 // densityDiffs[0] = 0.0001; //air
567
568 //--- Calculate the average material density for each material/density bin
569 std::map< std::pair<G4Material*,G4int>, matInfo* > newMateDens;
570
571 //---- Read the material densities
572 G4double dens;
573 for( G4int iz = 0; iz < fNoVoxelsZ; ++iz ) {
574 for( G4int iy = 0; iy < fNoVoxelsY; ++iy ) {
575 for( G4int ix = 0; ix < fNoVoxelsX; ++ix ) {
576 fin >> dens;
577 G4int copyNo = ix + (iy)*fNoVoxelsX + (iz)*fNoVoxelsX*fNoVoxelsY;
578
579 if( densityDiff != -1. ) continue;
580
581 //--- store the minimum and maximum density for each material
582 mpite = densiMinMax.find( G4int(fMateIDs[copyNo]) );
583 if( dens < (*mpite).second.first ) (*mpite).second.first = dens;
584 if( dens > (*mpite).second.second ) (*mpite).second.second = dens;
585 //--- Get material from original list of material in file
586 G4int mateID = G4int(fMateIDs[copyNo]);
587 std::map<G4int,G4Material*>::const_iterator imite =
588 fPhantomMaterialsOriginal.find(mateID);
589
590 //--- Check if density is equal to the original material density
591 if(std::fabs(dens - (*imite).second->GetDensity()/CLHEP::g*CLHEP::cm3 )
592 < 1.e-9 ) continue;
593
594 //--- Build material name with fPhantomMaterialsOriginal name+density
595 G4int densityBin = (G4int(dens/densityDiffs[mateID]));
596
597 G4String mateName = (*imite).second->GetName()
598 + G4UIcommand::ConvertToString(densityBin);
599 //--- Look if it is the first voxel with this material/densityBin
600 std::pair<G4Material*,G4int> matdens((*imite).second, densityBin );
601
602 auto mppite = newMateDens.find( matdens );
603 if( mppite != newMateDens.cend() ){
604 matInfo* mi = (*mppite).second;
605 mi->fSumdens += dens;
606 mi->fNvoxels++;
607 fMateIDs[copyNo] = fPhantomMaterialsOriginal.size()-1 + mi->fId;
608 } else {
609 matInfo* mi = new matInfo;
610 mi->fSumdens = dens;
611 mi->fNvoxels = 1;
612 mi->fId = G4int(newMateDens.size()+1);
613 newMateDens[matdens] = mi;
614 fMateIDs[copyNo] = fPhantomMaterialsOriginal.size()-1 + mi->fId;
615 }
616 }
617 }
618 }
619
620 if( densityDiff != -1. ) {
621 for( mpite = densiMinMax.begin(); mpite != densiMinMax.end(); ++mpite )
622 {
623#ifdef G4VERBOSE
624 G4cout << "DicomDetectorConstruction::ReadVoxelDensities"
625 << " ORIG MATERIALS DENSITY "
626 << (*mpite).first << " MIN " << (*mpite).second.first << " MAX "
627 << (*mpite).second.second << G4endl;
628#endif
629 }
630 }
631
632 //----- Build the list of phantom materials that go to Parameterisation
633 //--- Add original materials
634 for( auto mimite = fPhantomMaterialsOriginal.cbegin();
635 mimite != fPhantomMaterialsOriginal.cend(); ++mimite )
636 {
637 fMaterials.push_back( (*mimite).second );
638 }
639 //
640 //---- Build and add new materials
641 for( auto mppite= newMateDens.cbegin(); mppite!=newMateDens.cend(); ++mppite )
642 {
643 G4double averdens = (*mppite).second->fSumdens/(*mppite).second->fNvoxels;
644 G4double saverdens = G4int(1000.001*averdens)/1000.;
645#ifndef G4VERBOSE
646 G4cout << "DicomDetectorConstruction::ReadVoxelDensities AVER DENS "
647 << averdens << " -> "
648 << saverdens << " -> " << G4int(1000*averdens) << " "
649 << G4int(1000*averdens)/1000
650 << " " << G4int(1000*averdens)/1000. << G4endl;
651#endif
652
653 G4String mateName = ((*mppite).first).first->GetName() + "_"
654 + G4UIcommand::ConvertToString(saverdens);
656 (*mppite).first.first, G4float(averdens), mateName ) );
657 }
658}
std::vector< G4Material * > fMaterials
G4Material * BuildMaterialWithChangingDensity(const G4Material *origMate, G4float density, G4String newMateName)
Dicom detector construction.

◆ ReadPhantomDataFile()

void DicomDetectorConstruction::ReadPhantomDataFile ( const G4String fname)
protected

Definition at line 695 of file DicomDetectorConstruction.cc.

696{
697 G4cout << " DicomDetectorConstruction::ReadPhantomDataFile opening file "
698 << fname << G4endl; //GDEB
699
700#ifdef G4VERBOSE
701 G4cout << " DicomDetectorConstruction::ReadPhantomDataFile opening file "
702 << fname << G4endl;
703#endif
704
705 std::ifstream fin(fname.c_str(), std::ios_base::in);
706 if( !fin.is_open() ) {
707 G4Exception("DicomDetectorConstruction::ReadPhantomDataFile",
708 "",
709 FatalErrorInArgument,
710 G4String("File not found " + fname ).c_str());
711 }
712 //----- Define density differences (maximum density difference to create
713 // a new material)
714 char* part = std::getenv( "DICOM_CHANGE_MATERIAL_DENSITY" );
715 G4double densityDiff = -1.;
716 if( part ) densityDiff = G4UIcommand::ConvertToDouble(part);
717 if( densityDiff != -1. )
718 {
719 for( unsigned int ii = 0; ii < fOriginalMaterials.size(); ++ii )
720 {
721 fDensityDiffs[ii] = densityDiff; //currently all materials with
722 // same difference
723 }
724 }
725 else
726 {
727 if( fMaterials.size() == 0 ) { // do it only for first slice
728 for( unsigned int ii = 0; ii < fOriginalMaterials.size(); ++ii )
729 {
730 fMaterials.push_back( fOriginalMaterials[ii] );
731 }
732 }
733 }
734
735 //----- Read data header
736 DicomPhantomZSliceHeader* sliceHeader = new DicomPhantomZSliceHeader( fin );
737 fZSliceHeaders.push_back( sliceHeader );
738
739 //----- Read material indices
740 G4int nVoxels = sliceHeader->GetNoVoxels();
741
742 //--- If first slice, initiliaze fMateIDs
743 if( fZSliceHeaders.size() == 1 ) {
744 //fMateIDs = new unsigned int[fNoFiles*nVoxels];
745 fMateIDs = new size_t[fNoFiles*nVoxels];
746
747 }
748
749 unsigned int mateID;
750 // number of voxels from previously read slices
751 G4int voxelCopyNo = G4int((fZSliceHeaders.size()-1)*nVoxels);
752 for( G4int ii = 0; ii < nVoxels; ++ii, voxelCopyNo++ )
753 {
754 fin >> mateID;
755 fMateIDs[voxelCopyNo] = mateID;
756 }
757
758 //----- Read material densities and build new materials if two voxels have
759 // same material but its density is in a different density interval
760 // (size of density intervals defined by densityDiff)
761 G4double density;
762 // number of voxels from previously read slices
763 voxelCopyNo = G4int((fZSliceHeaders.size()-1)*nVoxels);
764 for( G4int ii = 0; ii < nVoxels; ++ii, voxelCopyNo++ )
765 {
766 fin >> density;
767
768 //-- Get material from list of original materials
769 mateID = unsigned(fMateIDs[voxelCopyNo]);
770 G4Material* mateOrig = fOriginalMaterials[mateID];
771
772 //-- Get density bin: middle point of the bin in which the current
773 // density is included
774 G4String newMateName = mateOrig->GetName();
775 G4float densityBin = 0.;
776 if( densityDiff != -1.) {
777 densityBin = G4float(fDensityDiffs[mateID]) *
778 (G4int(density/fDensityDiffs[mateID])+0.5);
779 //-- Build the new material name
780 newMateName += G4UIcommand::ConvertToString(densityBin);
781 }
782
783 //-- Look if a material with this name is already created
784 // (because a previous voxel was already in this density bin)
785 unsigned int im;
786 for( im = 0; im < fMaterials.size(); ++im )
787 {
788 if( fMaterials[im]->GetName() == newMateName ) {
789 break;
790 }
791 }
792 //-- If material is already created use index of this material
793 if( im != fMaterials.size() ) {
794 fMateIDs[voxelCopyNo] = im;
795 //-- else, create the material
796 } else {
797 if( densityDiff != -1.) {
798 fMaterials.push_back( BuildMaterialWithChangingDensity( mateOrig,
799 densityBin, newMateName ) );
800 fMateIDs[voxelCopyNo] = fMaterials.size()-1;
801 } else {
802 G4cerr << " im " << im << " < " << fMaterials.size() << " name "
803 << newMateName << G4endl;
804 G4Exception("DicomDetectorConstruction::ReadPhantomDataFile",
805 "",
806 FatalErrorInArgument,
807 "Wrong index in material"); //it should never reach here
808 }
809 }
810 }
811}
std::vector< DicomPhantomZSliceHeader * > fZSliceHeaders
std::map< G4int, G4double > fDensityDiffs
DicomPhantomZSliceHeader class.

◆ MergeZSliceHeaders()

void DicomDetectorConstruction::MergeZSliceHeaders ( )
protected

Definition at line 814 of file DicomDetectorConstruction.cc.

815{
816 //----- Images must have the same dimension ...
818 for( unsigned int ii = 1; ii < fZSliceHeaders.size(); ++ii )
819 {
821 }
822}

◆ BuildMaterialWithChangingDensity()

G4Material * DicomDetectorConstruction::BuildMaterialWithChangingDensity ( const G4Material origMate,
G4float  density,
G4String  newMateName 
)
protected

Definition at line 825 of file DicomDetectorConstruction.cc.

827{
828 //----- Copy original material, but with new density
829 G4int nelem = G4int(origMate->GetNumberOfElements());
830 G4Material* mate = new G4Material( newMateName, density*g/cm3, nelem,
831 kStateUndefined, STP_Temperature );
832
833 for( G4int ii = 0; ii < nelem; ++ii )
834 {
835 G4double frac = origMate->GetFractionVector()[ii];
836 G4Element* elem = const_cast<G4Element*>(origMate->GetElement(ii));
837 mate->AddElement( elem, frac );
838 }
839
840 return mate;
841}

◆ ConstructPhantomContainer()

void DicomDetectorConstruction::ConstructPhantomContainer ( )
protected

Definition at line 844 of file DicomDetectorConstruction.cc.

845{
846 //---- Extract number of voxels and voxel dimensions
850
854#ifdef G4VERBOSE
855 G4cout << " fNoVoxelsX " << fNoVoxelsX << " fVoxelHalfDimX " << fVoxelHalfDimX
856 <<G4endl;
857 G4cout << " fNoVoxelsY " << fNoVoxelsY << " fVoxelHalfDimY " << fVoxelHalfDimY
858 <<G4endl;
859 G4cout << " fNoVoxelsZ " << fNoVoxelsZ << " fVoxelHalfDimZ " << fVoxelHalfDimZ
860 <<G4endl;
861 G4cout << " totalPixels " << fNoVoxelsX*fNoVoxelsY*fNoVoxelsZ << G4endl;
862#endif
863
864 //----- Define the volume that contains all the voxels
865 fContainer_solid = new G4Box("phantomContainer",fNoVoxelsX*fVoxelHalfDimX,
870 //the material is not important, it will be fully filled by the voxels
871 fMaterials[0],
872 "phantomContainer",
873 0, 0, 0 );
874 //--- Place it on the world
875 G4double fOffsetX = (fZSliceHeaderMerged->GetMaxX() +
877 G4double fOffsetY = (fZSliceHeaderMerged->GetMaxY() +
879 G4double fOffsetZ = (fZSliceHeaderMerged->GetMaxZ() +
881 G4ThreeVector posCentreVoxels(fOffsetX,fOffsetY,fOffsetZ);
882#ifdef G4VERBOSE
883 G4cout << " placing voxel container volume at " << posCentreVoxels << G4endl;
884#endif
886 new G4PVPlacement(0, // rotation
887 posCentreVoxels,
888 fContainer_logic, // The logic volume
889 "phantomContainer", // Name
890 fWorld_logic, // Mother
891 false, // No op. bool.
892 1); // Copy number
893}

◆ ConstructPhantomContainerNew()

void DicomDetectorConstruction::ConstructPhantomContainerNew ( )
protected

Definition at line 896 of file DicomDetectorConstruction.cc.

897{
898#ifdef G4_DCMTK
899 //---- Extract number of voxels and voxel dimensions
900#ifdef G4VERBOSE
901 G4cout << " fNoVoxelsX " << fNoVoxelsX << " fVoxelHalfDimX " << fVoxelHalfDimX
902 <<G4endl;
903 G4cout << " fNoVoxelsY " << fNoVoxelsY << " fVoxelHalfDimY " << fVoxelHalfDimY
904 <<G4endl;
905 G4cout << " fNoVoxelsZ " << fNoVoxelsZ << " fVoxelHalfDimZ " << fVoxelHalfDimZ
906 <<G4endl;
907 G4cout << " totalPixels " << fNoVoxelsX*fNoVoxelsY*fNoVoxelsZ << G4endl;
908#endif
909
910 //----- Define the volume that contains all the voxels
911 fContainer_solid = new G4Box("phantomContainer",fNoVoxelsX*fVoxelHalfDimX,
916 //the material is not important, it will be fully filled by the voxels
917 fMaterials[0],
918 "phantomContainer",
919 0, 0, 0 );
920
921 G4ThreeVector posCentreVoxels((fMinX+fMaxX)/2.,
922 (fMinY+fMaxY)/2.,
923 (fMinZ+fMaxZ)/2.);
924#ifdef G4VERBOSE
925 G4cout << " placing voxel container volume at " << posCentreVoxels << G4endl;
926#endif
928 new G4PVPlacement(0, // rotation
929 posCentreVoxels,
930 fContainer_logic, // The logic volume
931 "phantomContainer", // Name
932 fWorld_logic, // Mother
933 false, // No op. bool.
934 1); // Copy number
935#endif
936}

◆ ConstructPhantom()

virtual void DicomDetectorConstruction::ConstructPhantom ( )
protectedpure virtual

◆ SetScorer()

void DicomDetectorConstruction::SetScorer ( G4LogicalVolume voxel_logic)
protected

Definition at line 944 of file DicomDetectorConstruction.cc.

945{
946
947#ifdef G4VERBOSE
948 G4cout << "\t SET SCORER : " << voxel_logic->GetName() << G4endl;
949#endif
950
951 fScorers.insert(voxel_logic);
952}
std::set< G4LogicalVolume * > fScorers

◆ ConstructSDandField()

void DicomDetectorConstruction::ConstructSDandField ( )
protectedvirtual

Definition at line 956 of file DicomDetectorConstruction.cc.

957{
958
959#ifdef G4VERBOSE
960 G4cout << "\t CONSTRUCT SD AND FIELD" << G4endl;
961#endif
962
963 //G4SDManager* SDman = G4SDManager::GetSDMpointer();
964
965 //SDman->SetVerboseLevel(1);
966
967 //
968 // Sensitive Detector Name
969 G4String concreteSDname = "phantomSD";
970 std::vector<G4String> scorer_names;
971 scorer_names.push_back(concreteSDname);
972 //------------------------
973 // MultiFunctionalDetector
974 //------------------------
975 //
976 // Define MultiFunctionalDetector with name.
977 // declare MFDet as a MultiFunctionalDetector scorer
979 new G4MultiFunctionalDetector(concreteSDname);
980 G4SDManager::GetSDMpointer()->AddNewDetector( MFDet );
981 char* nest = std::getenv( "DICOM_NESTED_PARAM" );
982 if( nest && G4String(nest) == "1" ) {
983 G4VPrimitiveScorer* dosedep =
984 new G4PSDoseDeposit3D("DoseDeposit", fNoVoxelsZ, fNoVoxelsY, fNoVoxelsX, 0, 2, 1); // nested param replica indexing
985 // - the last 3 arguments correspond to the replica depth for Z, Y and X respectively
986 MFDet->RegisterPrimitive(dosedep);
987 } else {
988 G4VPrimitiveScorer* dosedep = new G4PSDoseDeposit("DoseDeposit"); // regular geometry
989 MFDet->RegisterPrimitive(dosedep);
990 }
991
992 for(auto ite = fScorers.cbegin(); ite != fScorers.cend(); ++ite)
993 {
994 SetSensitiveDetector(*ite, MFDet);
995 }
996}

Member Data Documentation

◆ fAir

G4Material* DicomDetectorConstruction::fAir
protected

Definition at line 112 of file DicomDetectorConstruction.hh.

◆ fWorld_solid

G4Box* DicomDetectorConstruction::fWorld_solid
protected

Definition at line 115 of file DicomDetectorConstruction.hh.

◆ fWorld_logic

G4LogicalVolume* DicomDetectorConstruction::fWorld_logic
protected

Definition at line 116 of file DicomDetectorConstruction.hh.

◆ fWorld_phys

G4VPhysicalVolume* DicomDetectorConstruction::fWorld_phys
protected

Definition at line 117 of file DicomDetectorConstruction.hh.

◆ fContainer_solid

G4Box* DicomDetectorConstruction::fContainer_solid
protected

Definition at line 119 of file DicomDetectorConstruction.hh.

◆ fContainer_logic

G4LogicalVolume* DicomDetectorConstruction::fContainer_logic
protected

Definition at line 120 of file DicomDetectorConstruction.hh.

◆ fContainer_phys

G4VPhysicalVolume* DicomDetectorConstruction::fContainer_phys
protected

Definition at line 121 of file DicomDetectorConstruction.hh.

◆ fNoFiles

G4int DicomDetectorConstruction::fNoFiles
protected

Definition at line 123 of file DicomDetectorConstruction.hh.

◆ fOriginalMaterials

std::vector<G4Material*> DicomDetectorConstruction::fOriginalMaterials
protected

Definition at line 124 of file DicomDetectorConstruction.hh.

◆ fMaterials

std::vector<G4Material*> DicomDetectorConstruction::fMaterials
protected

Definition at line 125 of file DicomDetectorConstruction.hh.

◆ fMateIDs

size_t* DicomDetectorConstruction::fMateIDs
protected

Definition at line 128 of file DicomDetectorConstruction.hh.

◆ fDensityDiffs

std::map<G4int,G4double> DicomDetectorConstruction::fDensityDiffs
protected

Definition at line 131 of file DicomDetectorConstruction.hh.

◆ fZSliceHeaders

std::vector<DicomPhantomZSliceHeader*> DicomDetectorConstruction::fZSliceHeaders
protected

Definition at line 135 of file DicomDetectorConstruction.hh.

◆ fZSliceHeaderMerged

DicomPhantomZSliceHeader* DicomDetectorConstruction::fZSliceHeaderMerged
protected

Definition at line 137 of file DicomDetectorConstruction.hh.

◆ fNoVoxelsX

G4int DicomDetectorConstruction::fNoVoxelsX
protected

Definition at line 140 of file DicomDetectorConstruction.hh.

◆ fNoVoxelsY

G4int DicomDetectorConstruction::fNoVoxelsY
protected

Definition at line 140 of file DicomDetectorConstruction.hh.

◆ fNoVoxelsZ

G4int DicomDetectorConstruction::fNoVoxelsZ
protected

Definition at line 140 of file DicomDetectorConstruction.hh.

◆ fVoxelHalfDimX

G4double DicomDetectorConstruction::fVoxelHalfDimX
protected

Definition at line 141 of file DicomDetectorConstruction.hh.

◆ fVoxelHalfDimY

G4double DicomDetectorConstruction::fVoxelHalfDimY
protected

Definition at line 141 of file DicomDetectorConstruction.hh.

◆ fVoxelHalfDimZ

G4double DicomDetectorConstruction::fVoxelHalfDimZ
protected

Definition at line 141 of file DicomDetectorConstruction.hh.

◆ fMinX

G4double DicomDetectorConstruction::fMinX
protected

Definition at line 142 of file DicomDetectorConstruction.hh.

◆ fMinY

G4double DicomDetectorConstruction::fMinY
protected

Definition at line 142 of file DicomDetectorConstruction.hh.

◆ fMinZ

G4double DicomDetectorConstruction::fMinZ
protected

Definition at line 142 of file DicomDetectorConstruction.hh.

◆ fMaxX

G4double DicomDetectorConstruction::fMaxX
protected

Definition at line 143 of file DicomDetectorConstruction.hh.

◆ fMaxY

G4double DicomDetectorConstruction::fMaxY
protected

Definition at line 143 of file DicomDetectorConstruction.hh.

◆ fMaxZ

G4double DicomDetectorConstruction::fMaxZ
protected

Definition at line 143 of file DicomDetectorConstruction.hh.

◆ fPhantomMaterialsOriginal

std::map<G4int,G4Material*> DicomDetectorConstruction::fPhantomMaterialsOriginal
protected

Definition at line 145 of file DicomDetectorConstruction.hh.

◆ fMergedSlices

DicomPhantomZSliceMerged* DicomDetectorConstruction::fMergedSlices
protected

Definition at line 149 of file DicomDetectorConstruction.hh.

◆ fScorers

std::set<G4LogicalVolume*> DicomDetectorConstruction::fScorers
protected

Definition at line 151 of file DicomDetectorConstruction.hh.

◆ fConstructed

G4bool DicomDetectorConstruction::fConstructed
protected

Definition at line 153 of file DicomDetectorConstruction.hh.


The documentation for this class was generated from the following files:

Applications | User Support | Publications | Collaboration