ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DicomDetectorConstruction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DicomDetectorConstruction.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
29 //
30 
31 #include "globals.hh"
32 
33 #include "G4Box.hh"
34 #include "G4LogicalVolume.hh"
35 #include "G4VPhysicalVolume.hh"
36 #include "G4PVPlacement.hh"
37 #include "G4Material.hh"
38 #include "G4Element.hh"
39 #include "G4UIcommand.hh"
40 #include "G4PhysicalConstants.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 
48 #ifdef G4_DCMTK
49 #include "DicomFileMgr.hh"
50 #endif
51 #include "G4VisAttributes.hh"
52 
53 using CLHEP::m;
54 using CLHEP::cm3;
55 using CLHEP::mole;
56 using CLHEP::g;
57 using CLHEP::mg;
58 using CLHEP::perCent;
59 
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
61 DicomDetectorConstruction::DicomDetectorConstruction()
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  fNVoxelX(0),
79  fNVoxelY(0),
80  fNVoxelZ(0),
81  fVoxelHalfDimX(0),
82  fVoxelHalfDimY(0),
83  fVoxelHalfDimZ(0),
84 
85  fConstructed(false)
86 {
87 }
88 
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
90 DicomDetectorConstruction::~DicomDetectorConstruction()
91 {
92 }
93 
94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
95 G4VPhysicalVolume* DicomDetectorConstruction::Construct()
96 {
97  if(!fConstructed || fWorld_phys == 0) {
98  fConstructed = true;
99  InitialisationOfMaterials();
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 
111  fWorld_logic = new G4LogicalVolume( fWorld_solid,
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 );
123 
124 #ifdef G4_DCMTK
125  ReadPhantomDataNew();
126  ConstructPhantomContainerNew();
127 #else
128  ReadPhantomData();
129  ConstructPhantomContainer();
130 #endif
131 
132  ConstructPhantom();
133  }
134  return fWorld_phys;
135 }
136 
137 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........
138 void DicomDetectorConstruction::InitialisationOfMaterials()
139 {
140  // Creating elements :
141  G4double z, a, density;
142  G4String name, symbol;
143 
144  G4Element* elC = new G4Element( name = "Carbon",
145  symbol = "C",
146  z = 6.0, a = 12.011 * g/mole );
147  G4Element* elH = new G4Element( name = "Hydrogen",
148  symbol = "H",
149  z = 1.0, a = 1.008 * g/mole );
150  G4Element* elN = new G4Element( name = "Nitrogen",
151  symbol = "N",
152  z = 7.0, a = 14.007 * g/mole );
153  G4Element* elO = new G4Element( name = "Oxygen",
154  symbol = "O",
155  z = 8.0, a = 16.00 * g/mole );
156  G4Element* elNa = new G4Element( name = "Sodium",
157  symbol = "Na",
158  z= 11.0, a = 22.98977* g/mole );
159  G4Element* elMg = new G4Element( name = "Magnesium",
160  symbol = "Mg",
161  z = 12.0, a = 24.3050* g/mole );
162  G4Element* elP = new G4Element( name = "Phosphorus",
163  symbol = "P",
164  z = 15.0, a = 30.973976* g/mole );
165  G4Element* elS = new G4Element( name = "Sulfur",
166  symbol = "S",
167  z = 16.0,a = 32.065* g/mole );
168  G4Element* elCl = new G4Element( name = "Chlorine",
169  symbol = "P",
170  z = 17.0, a = 35.453* g/mole );
171  G4Element* elK = new G4Element( name = "Potassium",
172  symbol = "P",
173  z = 19.0, a = 30.0983* g/mole );
174 
175  G4Element* elFe = new G4Element( name = "Iron",
176  symbol = "Fe",
177  z = 26, a = 56.845* g/mole );
178 
179  G4Element* elCa = new G4Element( name="Calcium",
180  symbol = "Ca",
181  z = 20.0, a = 40.078* g/mole );
182 
183  G4Element* elZn = new G4Element( name = "Zinc",
184  symbol = "Zn",
185  z = 30.0,a = 65.382* g/mole );
186 
187  // Creating Materials :
188  G4int numberofElements;
189 
190  // Air
191  fAir = new G4Material( "Air",
192  1.290*mg/cm3,
193  numberofElements = 2 );
194  fAir->AddElement(elN, 0.7);
195  fAir->AddElement(elO, 0.3);
196 
197 
198  // Soft tissue (ICRP - NIST)
199  G4Material* softTissue = new G4Material ("SoftTissue", 1.00*g/cm3,
200  numberofElements = 13);
201  softTissue->AddElement(elH, 10.4472*perCent);
202  softTissue->AddElement(elC, 23.219*perCent);
203  softTissue->AddElement(elN, 2.488*perCent);
204  softTissue->AddElement(elO, 63.0238*perCent);
205  softTissue->AddElement(elNa, 0.113*perCent);
206  softTissue->AddElement(elMg, 0.0113*perCent);
207  softTissue->AddElement(elP, 0.113*perCent);
208  softTissue->AddElement(elS, 0.199*perCent);
209  softTissue->AddElement(elCl, 0.134*perCent);
210  softTissue->AddElement(elK, 0.199*perCent);
211  softTissue->AddElement(elCa, 0.023*perCent);
212  softTissue->AddElement(elFe, 0.005*perCent);
213  softTissue->AddElement(elZn, 0.003*perCent);
214 
215  // Lung Inhale
216  G4Material* lunginhale = new G4Material( "LungInhale",
217  density = 0.217*g/cm3,
218  numberofElements = 9);
219  lunginhale->AddElement(elH,0.103);
220  lunginhale->AddElement(elC,0.105);
221  lunginhale->AddElement(elN,0.031);
222  lunginhale->AddElement(elO,0.749);
223  lunginhale->AddElement(elNa,0.002);
224  lunginhale->AddElement(elP,0.002);
225  lunginhale->AddElement(elS,0.003);
226  lunginhale->AddElement(elCl,0.002);
227  lunginhale->AddElement(elK,0.003);
228 
229  // Lung exhale
230  G4Material* lungexhale = new G4Material( "LungExhale",
231  density = 0.508*g/cm3,
232  numberofElements = 9 );
233  lungexhale->AddElement(elH,0.103);
234  lungexhale->AddElement(elC,0.105);
235  lungexhale->AddElement(elN,0.031);
236  lungexhale->AddElement(elO,0.749);
237  lungexhale->AddElement(elNa,0.002);
238  lungexhale->AddElement(elP,0.002);
239  lungexhale->AddElement(elS,0.003);
240  lungexhale->AddElement(elCl,0.002);
241  lungexhale->AddElement(elK,0.003);
242 
243  // Adipose tissue
244  G4Material* adiposeTissue = new G4Material( "AdiposeTissue",
245  density = 0.967*g/cm3,
246  numberofElements = 7);
247  adiposeTissue->AddElement(elH,0.114);
248  adiposeTissue->AddElement(elC,0.598);
249  adiposeTissue->AddElement(elN,0.007);
250  adiposeTissue->AddElement(elO,0.278);
251  adiposeTissue->AddElement(elNa,0.001);
252  adiposeTissue->AddElement(elS,0.001);
253  adiposeTissue->AddElement(elCl,0.001);
254 
255  // Brain (ICRP - NIST)
256  G4Material* brainTissue = new G4Material ("BrainTissue", 1.03 * g/cm3,
257  numberofElements = 13);
258  brainTissue->AddElement(elH, 11.0667*perCent);
259  brainTissue->AddElement(elC, 12.542*perCent);
260  brainTissue->AddElement(elN, 1.328*perCent);
261  brainTissue->AddElement(elO, 73.7723*perCent);
262  brainTissue->AddElement(elNa, 0.1840*perCent);
263  brainTissue->AddElement(elMg, 0.015*perCent);
264  brainTissue->AddElement(elP, 0.356*perCent);
265  brainTissue->AddElement(elS, 0.177*perCent);
266  brainTissue->AddElement(elCl, 0.236*perCent);
267  brainTissue->AddElement(elK, 0.31*perCent);
268  brainTissue->AddElement(elCa, 0.009*perCent);
269  brainTissue->AddElement(elFe, 0.005*perCent);
270  brainTissue->AddElement(elZn, 0.001*perCent);
271 
272 
273  // Breast
274  G4Material* breast = new G4Material( "Breast",
275  density = 0.990*g/cm3,
276  numberofElements = 8 );
277  breast->AddElement(elH,0.109);
278  breast->AddElement(elC,0.506);
279  breast->AddElement(elN,0.023);
280  breast->AddElement(elO,0.358);
281  breast->AddElement(elNa,0.001);
282  breast->AddElement(elP,0.001);
283  breast->AddElement(elS,0.001);
284  breast->AddElement(elCl,0.001);
285 
286  // Spinal Disc
287  G4Material* spinalDisc = new G4Material ("SpinalDisc", 1.10 * g/cm3,
288  numberofElements = 8);
289  spinalDisc->AddElement(elH, 9.60*perCent);
290  spinalDisc->AddElement(elC, 9.90*perCent);
291  spinalDisc->AddElement(elN, 2.20*perCent);
292  spinalDisc->AddElement(elO, 74.40*perCent);
293  spinalDisc->AddElement(elNa, 0.50*perCent);
294  spinalDisc->AddElement(elP, 2.20*perCent);
295  spinalDisc->AddElement(elS, 0.90*perCent);
296  spinalDisc->AddElement(elCl, 0.30*perCent);
297 
298 
299  // Water
300  G4Material* water = new G4Material( "Water",
301  density = 1.0*g/cm3,
302  numberofElements = 2 );
303  water->AddElement(elH,0.112);
304  water->AddElement(elO,0.888);
305 
306  // Muscle
307  G4Material* muscle = new G4Material( "Muscle",
308  density = 1.061*g/cm3,
309  numberofElements = 9 );
310  muscle->AddElement(elH,0.102);
311  muscle->AddElement(elC,0.143);
312  muscle->AddElement(elN,0.034);
313  muscle->AddElement(elO,0.710);
314  muscle->AddElement(elNa,0.001);
315  muscle->AddElement(elP,0.002);
316  muscle->AddElement(elS,0.003);
317  muscle->AddElement(elCl,0.001);
318  muscle->AddElement(elK,0.004);
319 
320  // Liver
321  G4Material* liver = new G4Material( "Liver",
322  density = 1.071*g/cm3,
323  numberofElements = 9);
324  liver->AddElement(elH,0.102);
325  liver->AddElement(elC,0.139);
326  liver->AddElement(elN,0.030);
327  liver->AddElement(elO,0.716);
328  liver->AddElement(elNa,0.002);
329  liver->AddElement(elP,0.003);
330  liver->AddElement(elS,0.003);
331  liver->AddElement(elCl,0.002);
332  liver->AddElement(elK,0.003);
333 
334  // Tooth Dentin
335  G4Material* toothDentin = new G4Material ("ToothDentin", 2.14 * g/cm3,
336  numberofElements = 10);
337  toothDentin->AddElement(elH, 2.67*perCent);
338  toothDentin->AddElement(elC, 12.77*perCent);
339  toothDentin->AddElement(elN, 4.27*perCent);
340  toothDentin->AddElement(elO, 40.40*perCent);
341  toothDentin->AddElement(elNa, 0.65*perCent);
342  toothDentin->AddElement(elMg, 0.59*perCent);
343  toothDentin->AddElement(elP, 11.86*perCent);
344  toothDentin->AddElement(elCl, 0.04*perCent);
345  toothDentin->AddElement(elCa, 26.74*perCent);
346  toothDentin->AddElement(elZn, 0.01*perCent);
347 
348 
349  // Trabecular Bone
350  G4Material* trabecularBone = new G4Material("TrabecularBone",
351  density = 1.159*g/cm3,
352  numberofElements = 12 );
353  trabecularBone->AddElement(elH,0.085);
354  trabecularBone->AddElement(elC,0.404);
355  trabecularBone->AddElement(elN,0.058);
356  trabecularBone->AddElement(elO,0.367);
357  trabecularBone->AddElement(elNa,0.001);
358  trabecularBone->AddElement(elMg,0.001);
359  trabecularBone->AddElement(elP,0.034);
360  trabecularBone->AddElement(elS,0.002);
361  trabecularBone->AddElement(elCl,0.002);
362  trabecularBone->AddElement(elK,0.001);
363  trabecularBone->AddElement(elCa,0.044);
364  trabecularBone->AddElement(elFe,0.001);
365 
366  // Trabecular bone used in the DICOM Head
367 
368  G4Material* trabecularBone_head = new G4Material ("TrabecularBone_HEAD",
369  1.18 * g/cm3,
370  numberofElements = 12);
371  trabecularBone_head->AddElement(elH, 8.50*perCent);
372  trabecularBone_head->AddElement(elC, 40.40*perCent);
373  trabecularBone_head->AddElement(elN, 2.80*perCent);
374  trabecularBone_head->AddElement(elO, 36.70*perCent);
375  trabecularBone_head->AddElement(elNa, 0.10*perCent);
376  trabecularBone_head->AddElement(elMg, 0.10*perCent);
377  trabecularBone_head->AddElement(elP, 3.40*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);
427 
428 #ifdef DICOM_USE_HEAD
429  //----- Put the materials in a vector HEAD PHANTOM
430  fOriginalMaterials.push_back(fAir); //0.00129 g/cm3
431  fOriginalMaterials.push_back(softTissue); // 1.055 g/cm3
432  fOriginalMaterials.push_back(brainTissue); // 1.07 g/cm3
433  fOriginalMaterials.push_back(spinalDisc); // 1.10 g/cm3
434  fOriginalMaterials.push_back(trabecularBone_head); // 1.13 g/cm3
435  fOriginalMaterials.push_back(toothDentin); // 1.66 g/cm3
436  fOriginalMaterials.push_back(corticalBone); // 1.75 g/cm3
437  fOriginalMaterials.push_back(toothEnamel); // 2.04 g/cm3
438  G4cout << "The materials of the DICOM Head have been used" << G4endl;
439 #else
440  fOriginalMaterials.push_back(fAir); // rho = 0.00129
441  fOriginalMaterials.push_back(lunginhale); // rho = 0.217
442  fOriginalMaterials.push_back(lungexhale); // rho = 0.508
443  fOriginalMaterials.push_back(adiposeTissue); // rho = 0.967
444  fOriginalMaterials.push_back(breast ); // rho = 0.990
445  fOriginalMaterials.push_back(water); // rho = 1.018
446  fOriginalMaterials.push_back(muscle); // rho = 1.061
447  fOriginalMaterials.push_back(liver); // rho = 1.071
448  fOriginalMaterials.push_back(trabecularBone); // rho = 1.159 - HEAD PHANTOM
449  fOriginalMaterials.push_back(denseBone); // rho = 1.575
450  G4cout << "Default materials of the DICOM Extended examples have been used"
451  << G4endl;
452 #endif
453 }
454 
455 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
456 void DicomDetectorConstruction::ReadPhantomDataNew()
457 {
458 #ifdef G4_DCMTK
459  G4String fileName = DicomFileMgr::GetInstance()->GetFileOutName();
460 
461  std::ifstream fin(fileName);
462  std::vector<G4String> wl;
463  G4int nMaterials;
464  fin >> nMaterials;
465  G4String mateName;
466  G4int nmate;
467  for( G4int ii = 0; ii < nMaterials; ++ii )
468  {
469  fin >> nmate;
470  fin >> mateName;
471  if( mateName[0] == '"' && mateName[mateName.length()-1] == '"' ) {
472  mateName = mateName.substr(1,mateName.length()-2);
473  }
474  G4cout << "GmReadPhantomG4Geometry::ReadPhantomData reading nmate "
475  << ii << " = " << nmate
476  << " mate " << mateName << G4endl;
477  if( ii != nmate )
478  G4Exception("GmReadPhantomG4Geometry::ReadPhantomData",
479  "Wrong argument", FatalErrorInArgument,
480  "Material number should be in increasing order:wrong material number");
481 
482  G4Material* mate = 0;
483  const G4MaterialTable* matTab = G4Material::GetMaterialTable();
484  for( auto matite = matTab->cbegin(); matite != matTab->cend(); ++matite )
485  {
486  if( (*matite)->GetName() == mateName ) {
487  mate = *matite;
488  }
489  }
490  if( mate == 0 ) {
491  mate = G4NistManager::Instance()->FindOrBuildMaterial(mateName);
492  }
493  if( !mate ) G4Exception("GmReadPhantomG4Geometry::ReadPhantomData",
494  "Wrong argument",
495  FatalErrorInArgument,
496  ("Material not found" + mateName).c_str());
497  thePhantomMaterialsOriginal[nmate] = mate;
498  }
499 
500  fin >> fNVoxelX >> fNVoxelY >> fNVoxelZ;
501  G4cout << "GmReadPhantomG4Geometry::ReadPhantomData fNVoxel X/Y/Z "
502  << fNVoxelX << " "
503  << fNVoxelY << " " << fNVoxelZ << G4endl;
504  fin >> fMinX >> fMaxX;
505  fin >> fMinY >> fMaxY;
506  fin >> fMinZ >> fMaxZ;
507  fVoxelHalfDimX = (fMaxX-fMinX)/fNVoxelX/2.;
508  fVoxelHalfDimY = (fMaxY-fMinY)/fNVoxelY/2.;
509  fVoxelHalfDimZ = (fMaxZ-fMinZ)/fNVoxelZ/2.;
510 #ifdef G4VERBOSE
511  G4cout << " Extension in X " << fMinX << " " << fMaxX << G4endl
512  << " Extension in Y " << fMinY << " " << fMaxY << G4endl
513  << " Extension in Z " << fMinZ << " " << fMaxZ << G4endl;
514 #endif
515 
516  fMateIDs = new size_t[fNVoxelX*fNVoxelY*fNVoxelZ];
517  for( G4int iz = 0; iz < fNVoxelZ; ++iz ) {
518  for( G4int iy = 0; iy < fNVoxelY; ++iy ) {
519  for( G4int ix = 0; ix < fNVoxelX; ++ix ) {
520  G4int mateID;
521  fin >> mateID;
522  G4int nnew = ix + (iy)*fNVoxelX + (iz)*fNVoxelX*fNVoxelY;
523  if( mateID < 0 || mateID >= nMaterials ) {
524  G4Exception("GmReadPhantomG4Geometry::ReadPhantomData",
525  "Wrong index in phantom file",
526  FatalException,
527  G4String("It should be between 0 and "
528  + G4UIcommand::ConvertToString(nMaterials-1)
529  + ", while it is "
530  + G4UIcommand::ConvertToString(mateID)).c_str());
531  }
532  fMateIDs[nnew] = mateID;
533  }
534  }
535  }
536 
537  ReadVoxelDensities( fin );
538 
539  fin.close();
540 #endif
541 
542 }
543 
544 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
545 void DicomDetectorConstruction::ReadVoxelDensities( std::ifstream& fin )
546 {
547  G4String stemp;
548  std::map<G4int, std::pair<G4double,G4double> > densiMinMax;
549  std::map<G4int, std::pair<G4double,G4double> >::iterator mpite;
550  for( G4int ii = 0; ii < G4int(thePhantomMaterialsOriginal.size()); ++ii )
551  {
552  densiMinMax[ii] = std::pair<G4double,G4double>(DBL_MAX,-DBL_MAX);
553  }
554 
555  char* part = std::getenv( "DICOM_CHANGE_MATERIAL_DENSITY" );
556  G4double densityDiff = -1.;
557  if( part ) densityDiff = G4UIcommand::ConvertToDouble(part);
558 
559  std::map<G4int,G4double> densityDiffs;
560  for( G4int ii = 0; ii < G4int(thePhantomMaterialsOriginal.size()); ++ii )
561  {
562  densityDiffs[ii] = densityDiff; //currently all materials with same step
563  }
564  // densityDiffs[0] = 0.0001; //air
565 
566  //--- Calculate the average material density for each material/density bin
567  std::map< std::pair<G4Material*,G4int>, matInfo* > newMateDens;
568 
569  //---- Read the material densities
570  G4double dens;
571  for( G4int iz = 0; iz < fNVoxelZ; ++iz ) {
572  for( G4int iy = 0; iy < fNVoxelY; ++iy ) {
573  for( G4int ix = 0; ix < fNVoxelX; ++ix ) {
574  fin >> dens;
575  G4int copyNo = ix + (iy)*fNVoxelX + (iz)*fNVoxelX*fNVoxelY;
576 
577  if( densityDiff != -1. ) continue;
578 
579  //--- store the minimum and maximum density for each material
580  mpite = densiMinMax.find( G4int(fMateIDs[copyNo]) );
581  if( dens < (*mpite).second.first ) (*mpite).second.first = dens;
582  if( dens > (*mpite).second.second ) (*mpite).second.second = dens;
583  //--- Get material from original list of material in file
584  G4int mateID = G4int(fMateIDs[copyNo]);
585  std::map<G4int,G4Material*>::const_iterator imite =
586  thePhantomMaterialsOriginal.find(mateID);
587 
588  //--- Check if density is equal to the original material density
589  if(std::fabs(dens - (*imite).second->GetDensity()/CLHEP::g*CLHEP::cm3 )
590  < 1.e-9 ) continue;
591 
592  //--- Build material name with thePhantomMaterialsOriginal name+density
593  G4int densityBin = (G4int(dens/densityDiffs[mateID]));
594 
595  G4String mateName = (*imite).second->GetName()
596  + G4UIcommand::ConvertToString(densityBin);
597  //--- Look if it is the first voxel with this material/densityBin
598  std::pair<G4Material*,G4int> matdens((*imite).second, densityBin );
599 
600  auto mppite = newMateDens.find( matdens );
601  if( mppite != newMateDens.cend() ){
602  matInfo* mi = (*mppite).second;
603  mi->fSumdens += dens;
604  mi->fNvoxels++;
605  fMateIDs[copyNo] = thePhantomMaterialsOriginal.size()-1 + mi->fId;
606  } else {
607  matInfo* mi = new matInfo;
608  mi->fSumdens = dens;
609  mi->fNvoxels = 1;
610  mi->fId = G4int(newMateDens.size()+1);
611  newMateDens[matdens] = mi;
612  fMateIDs[copyNo] = thePhantomMaterialsOriginal.size()-1 + mi->fId;
613  }
614  }
615  }
616  }
617 
618  if( densityDiff != -1. ) {
619  for( mpite = densiMinMax.begin(); mpite != densiMinMax.end(); ++mpite )
620  {
621 #ifdef G4VERBOSE
622  G4cout << "DicomDetectorConstruction::ReadVoxelDensities"
623  << " ORIG MATERIALS DENSITY "
624  << (*mpite).first << " MIN " << (*mpite).second.first << " MAX "
625  << (*mpite).second.second << G4endl;
626 #endif
627  }
628  }
629 
630  //----- Build the list of phantom materials that go to Parameterisation
631  //--- Add original materials
632  for( auto mimite = thePhantomMaterialsOriginal.cbegin();
633  mimite != thePhantomMaterialsOriginal.cend(); ++mimite )
634  {
635  fMaterials.push_back( (*mimite).second );
636  }
637  //
638  //---- Build and add new materials
639  for( auto mppite= newMateDens.cbegin(); mppite!=newMateDens.cend(); ++mppite )
640  {
641  G4double averdens = (*mppite).second->fSumdens/(*mppite).second->fNvoxels;
642  G4double saverdens = G4int(1000.001*averdens)/1000.;
643 #ifndef G4VERBOSE
644  G4cout << "DicomDetectorConstruction::ReadVoxelDensities AVER DENS "
645  << averdens << " -> "
646  << saverdens << " -> " << G4int(1000*averdens) << " "
647  << G4int(1000*averdens)/1000
648  << " " << G4int(1000*averdens)/1000. << G4endl;
649 #endif
650 
651  G4String mateName = ((*mppite).first).first->GetName() + "_"
652  + G4UIcommand::ConvertToString(saverdens);
653  fMaterials.push_back( BuildMaterialWithChangingDensity(
654  (*mppite).first.first, G4float(averdens), mateName ) );
655  }
656 }
657 
658 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.......
659 void DicomDetectorConstruction::ReadPhantomData()
660 {
661  G4String dataFile = DicomHandler::GetDicomDataFile();
662  std::ifstream finDF(dataFile.c_str());
663  G4String fname;
664 
665 if(finDF.good() != 1 )
666  {
667  G4String descript = "Problem reading data file: "+dataFile;
668  G4Exception(" DicomDetectorConstruction::ReadPhantomData"," ",
669  FatalException,descript);
670  }
671 
672  G4int compression;
673  finDF >> compression; // not used here
674  finDF >> fNoFiles;
675 
676  for(G4int i = 0; i < fNoFiles; ++i )
677  {
678 
679  finDF >> fname;
680 
681  //--- Read one data file
682  fname += ".g4dcm";
683 
684  ReadPhantomDataFile(fname);
685  }
686 
687  //----- Merge data headers
688  MergeZSliceHeaders();
689  finDF.close();
690 }
691 
692 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........
693 void DicomDetectorConstruction::ReadPhantomDataFile(const G4String& fname)
694 {
695  G4cout << " DicomDetectorConstruction::ReadPhantomDataFile opening file "
696  << fname << G4endl; //GDEB
697 
698 #ifdef G4VERBOSE
699  G4cout << " DicomDetectorConstruction::ReadPhantomDataFile opening file "
700  << fname << G4endl;
701 #endif
702 
703  std::ifstream fin(fname.c_str(), std::ios_base::in);
704  if( !fin.is_open() ) {
705  G4Exception("DicomDetectorConstruction::ReadPhantomDataFile",
706  "",
707  FatalErrorInArgument,
708  G4String("File not found " + fname ).c_str());
709  }
710  //----- Define density differences (maximum density difference to create
711  // a new material)
712  char* part = std::getenv( "DICOM_CHANGE_MATERIAL_DENSITY" );
713  G4double densityDiff = -1.;
714  if( part ) densityDiff = G4UIcommand::ConvertToDouble(part);
715  if( densityDiff != -1. )
716  {
717  for( unsigned int ii = 0; ii < fOriginalMaterials.size(); ++ii )
718  {
719  fDensityDiffs[ii] = densityDiff; //currently all materials with
720  // same difference
721  }
722  }
723  else
724  {
725  if( fMaterials.size() == 0 ) { // do it only for first slice
726  for( unsigned int ii = 0; ii < fOriginalMaterials.size(); ++ii )
727  {
728  fMaterials.push_back( fOriginalMaterials[ii] );
729  }
730  }
731  }
732 
733  //----- Read data header
734  DicomPhantomZSliceHeader* sliceHeader = new DicomPhantomZSliceHeader( fin );
735  fZSliceHeaders.push_back( sliceHeader );
736 
737  //----- Read material indices
738  G4int nVoxels = sliceHeader->GetNoVoxels();
739 
740  //--- If first slice, initiliaze fMateIDs
741  if( fZSliceHeaders.size() == 1 ) {
742  //fMateIDs = new unsigned int[fNoFiles*nVoxels];
743  fMateIDs = new size_t[fNoFiles*nVoxels];
744 
745  }
746 
747  unsigned int mateID;
748  // number of voxels from previously read slices
749  G4int voxelCopyNo = G4int((fZSliceHeaders.size()-1)*nVoxels);
750  for( G4int ii = 0; ii < nVoxels; ++ii, voxelCopyNo++ )
751  {
752  fin >> mateID;
753  fMateIDs[voxelCopyNo] = mateID;
754  }
755 
756  //----- Read material densities and build new materials if two voxels have
757  // same material but its density is in a different density interval
758  // (size of density intervals defined by densityDiff)
759  G4double density;
760  // number of voxels from previously read slices
761  voxelCopyNo = G4int((fZSliceHeaders.size()-1)*nVoxels);
762  for( G4int ii = 0; ii < nVoxels; ++ii, voxelCopyNo++ )
763  {
764  fin >> density;
765 
766  //-- Get material from list of original materials
767  mateID = unsigned(fMateIDs[voxelCopyNo]);
768  G4Material* mateOrig = fOriginalMaterials[mateID];
769 
770  //-- Get density bin: middle point of the bin in which the current
771  // density is included
772  G4String newMateName = mateOrig->GetName();
773  G4float densityBin = 0.;
774  if( densityDiff != -1.) {
775  densityBin = G4float(fDensityDiffs[mateID]) *
776  (G4int(density/fDensityDiffs[mateID])+0.5);
777  //-- Build the new material name
778  newMateName += G4UIcommand::ConvertToString(densityBin);
779  }
780 
781  //-- Look if a material with this name is already created
782  // (because a previous voxel was already in this density bin)
783  unsigned int im;
784  for( im = 0; im < fMaterials.size(); ++im )
785  {
786  if( fMaterials[im]->GetName() == newMateName ) {
787  break;
788  }
789  }
790  //-- If material is already created use index of this material
791  if( im != fMaterials.size() ) {
792  fMateIDs[voxelCopyNo] = im;
793  //-- else, create the material
794  } else {
795  if( densityDiff != -1.) {
796  fMaterials.push_back( BuildMaterialWithChangingDensity( mateOrig,
797  densityBin, newMateName ) );
798  fMateIDs[voxelCopyNo] = fMaterials.size()-1;
799  } else {
800  G4cerr << " im " << im << " < " << fMaterials.size() << " name "
801  << newMateName << G4endl;
802  G4Exception("DicomDetectorConstruction::ReadPhantomDataFile",
803  "",
804  FatalErrorInArgument,
805  "Wrong index in material"); //it should never reach here
806  }
807  }
808  }
809 }
810 
811 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
812 void DicomDetectorConstruction::MergeZSliceHeaders()
813 {
814  //----- Images must have the same dimension ...
815  fZSliceHeaderMerged = new DicomPhantomZSliceHeader( *fZSliceHeaders[0] );
816  for( unsigned int ii = 1; ii < fZSliceHeaders.size(); ++ii )
817  {
818  *fZSliceHeaderMerged += *fZSliceHeaders[ii];
819  }
820 }
821 
822 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
823 G4Material* DicomDetectorConstruction::BuildMaterialWithChangingDensity(
824  const G4Material* origMate, G4float density, G4String newMateName )
825 {
826  //----- Copy original material, but with new density
827  G4int nelem = G4int(origMate->GetNumberOfElements());
828  G4Material* mate = new G4Material( newMateName, density*g/cm3, nelem,
829  kStateUndefined, STP_Temperature );
830 
831  for( G4int ii = 0; ii < nelem; ++ii )
832  {
833  G4double frac = origMate->GetFractionVector()[ii];
834  G4Element* elem = const_cast<G4Element*>(origMate->GetElement(ii));
835  mate->AddElement( elem, frac );
836  }
837 
838  return mate;
839 }
840 
841 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
842 void DicomDetectorConstruction::ConstructPhantomContainer()
843 {
844  //---- Extract number of voxels and voxel dimensions
845  fNVoxelX = fZSliceHeaderMerged->GetNoVoxelX();
846  fNVoxelY = fZSliceHeaderMerged->GetNoVoxelY();
847  fNVoxelZ = fZSliceHeaderMerged->GetNoVoxelZ();
848 
849  fVoxelHalfDimX = fZSliceHeaderMerged->GetVoxelHalfX();
850  fVoxelHalfDimY = fZSliceHeaderMerged->GetVoxelHalfY();
851  fVoxelHalfDimZ = fZSliceHeaderMerged->GetVoxelHalfZ();
852 #ifdef G4VERBOSE
853  G4cout << " fNVoxelX " << fNVoxelX << " fVoxelHalfDimX " << fVoxelHalfDimX
854  <<G4endl;
855  G4cout << " fNVoxelY " << fNVoxelY << " fVoxelHalfDimY " << fVoxelHalfDimY
856  <<G4endl;
857  G4cout << " fNVoxelZ " << fNVoxelZ << " fVoxelHalfDimZ " << fVoxelHalfDimZ
858  <<G4endl;
859  G4cout << " totalPixels " << fNVoxelX*fNVoxelY*fNVoxelZ << G4endl;
860 #endif
861 
862  //----- Define the volume that contains all the voxels
863  fContainer_solid = new G4Box("phantomContainer",fNVoxelX*fVoxelHalfDimX,
864  fNVoxelY*fVoxelHalfDimY,
865  fNVoxelZ*fVoxelHalfDimZ);
866  fContainer_logic =
867  new G4LogicalVolume( fContainer_solid,
868  //the material is not important, it will be fully filled by the voxels
869  fMaterials[0],
870  "phantomContainer",
871  0, 0, 0 );
872  //--- Place it on the world
873  G4double fOffsetX = (fZSliceHeaderMerged->GetMaxX() +
874  fZSliceHeaderMerged->GetMinX() ) /2.;
875  G4double fOffsetY = (fZSliceHeaderMerged->GetMaxY() +
876  fZSliceHeaderMerged->GetMinY() ) /2.;
877  G4double fOffsetZ = (fZSliceHeaderMerged->GetMaxZ() +
878  fZSliceHeaderMerged->GetMinZ() ) /2.;
879  G4ThreeVector posCentreVoxels(fOffsetX,fOffsetY,fOffsetZ);
880 #ifdef G4VERBOSE
881  G4cout << " placing voxel container volume at " << posCentreVoxels << G4endl;
882 #endif
883  fContainer_phys =
884  new G4PVPlacement(0, // rotation
885  posCentreVoxels,
886  fContainer_logic, // The logic volume
887  "phantomContainer", // Name
888  fWorld_logic, // Mother
889  false, // No op. bool.
890  1); // Copy number
891 }
892 
893 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
894 void DicomDetectorConstruction::ConstructPhantomContainerNew()
895 {
896 #ifdef G4_DCMTK
897  //---- Extract number of voxels and voxel dimensions
898 #ifdef G4VERBOSE
899  G4cout << " fNVoxelX " << fNVoxelX << " fVoxelHalfDimX " << fVoxelHalfDimX
900  <<G4endl;
901  G4cout << " fNVoxelY " << fNVoxelY << " fVoxelHalfDimY " << fVoxelHalfDimY
902  <<G4endl;
903  G4cout << " fNVoxelZ " << fNVoxelZ << " fVoxelHalfDimZ " << fVoxelHalfDimZ
904  <<G4endl;
905  G4cout << " totalPixels " << fNVoxelX*fNVoxelY*fNVoxelZ << G4endl;
906 #endif
907 
908  //----- Define the volume that contains all the voxels
909  fContainer_solid = new G4Box("phantomContainer",fNVoxelX*fVoxelHalfDimX,
910  fNVoxelY*fVoxelHalfDimY,
911  fNVoxelZ*fVoxelHalfDimZ);
912  fContainer_logic =
913  new G4LogicalVolume( fContainer_solid,
914  //the material is not important, it will be fully filled by the voxels
915  fMaterials[0],
916  "phantomContainer",
917  0, 0, 0 );
918 
919  G4ThreeVector posCentreVoxels((fMinX+fMaxX)/2.,
920  (fMinY+fMaxY)/2.,
921  (fMinZ+fMaxZ)/2.);
922 #ifdef G4VERBOSE
923  G4cout << " placing voxel container volume at " << posCentreVoxels << G4endl;
924 #endif
925  fContainer_phys =
926  new G4PVPlacement(0, // rotation
927  posCentreVoxels,
928  fContainer_logic, // The logic volume
929  "phantomContainer", // Name
930  fWorld_logic, // Mother
931  false, // No op. bool.
932  1); // Copy number
933 #endif
934 }
935 
936 #include "G4SDManager.hh"
937 #include "G4MultiFunctionalDetector.hh"
938 #include "G4PSDoseDeposit.hh"
939 #include "G4PSDoseDeposit3D.hh"
940 
941 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.
942 void DicomDetectorConstruction::SetScorer(G4LogicalVolume* voxel_logic)
943 {
944 
945 #ifdef G4VERBOSE
946  G4cout << "\t SET SCORER : " << voxel_logic->GetName() << G4endl;
947 #endif
948 
949  fScorers.insert(voxel_logic);
950 }
951 
952 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
953 
954 void DicomDetectorConstruction::ConstructSDandField()
955 {
956 
957 #ifdef G4VERBOSE
958  G4cout << "\t CONSTRUCT SD AND FIELD" << G4endl;
959 #endif
960 
961  //G4SDManager* SDman = G4SDManager::GetSDMpointer();
962 
963  //SDman->SetVerboseLevel(1);
964 
965  //
966  // Sensitive Detector Name
967  G4String concreteSDname = "phantomSD";
968  std::vector<G4String> scorer_names;
969  scorer_names.push_back(concreteSDname);
970  //------------------------
971  // MultiFunctionalDetector
972  //------------------------
973  //
974  // Define MultiFunctionalDetector with name.
975  // declare MFDet as a MultiFunctionalDetector scorer
976  G4MultiFunctionalDetector* MFDet =
977  new G4MultiFunctionalDetector(concreteSDname);
978  G4SDManager::GetSDMpointer()->AddNewDetector( MFDet );
979  //G4VPrimitiveScorer* dosedep = new G4PSDoseDeposit("DoseDeposit");
980  G4VPrimitiveScorer* dosedep =
981  new G4PSDoseDeposit3D("DoseDeposit", fNVoxelX, fNVoxelY, fNVoxelZ);
982  MFDet->RegisterPrimitive(dosedep);
983 
984  for(auto ite = fScorers.cbegin(); ite != fScorers.cend(); ++ite)
985  {
986  SetSensitiveDetector(*ite, MFDet);
987  }
988 }