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