Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DicomDetectorConstruction.cc
Go to the documentation of this file.
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 105190 2017-07-14 12:08:43Z 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"
44 #include "CLHEP/Units/SystemOfUnits.h"
45 
46 #include "DicomDetectorConstruction.hh"
47 #include "DicomPhantomZSliceHeader.hh"
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 
62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
63 DicomDetectorConstruction::DicomDetectorConstruction()
64  : G4VUserDetectorConstruction(),
65  fAir(0),
66 
67  fWorld_solid(0),
68  fWorld_logic(0),
69  fWorld_phys(0),
70 
71  fContainer_solid(0),
72  fContainer_logic(0),
73  fContainer_phys(0),
74 
75  fNoFiles(0),
76  fMateIDs(0),
77 
78  fZSliceHeaderMerged(0),
79 
80  fNVoxelX(0),
81  fNVoxelY(0),
82  fNVoxelZ(0),
83  fVoxelHalfDimX(0),
84  fVoxelHalfDimY(0),
85  fVoxelHalfDimZ(0),
86 
87  fConstructed(false)
88 {
89 
90 }
91 
92 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
93 DicomDetectorConstruction::~DicomDetectorConstruction()
94 {
95 }
96 
97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
98 G4VPhysicalVolume* DicomDetectorConstruction::Construct()
99 {
100  if(!fConstructed || fWorld_phys == 0) {
101  fConstructed = true;
102  InitialisationOfMaterials();
103 
104  //----- Build world
105  G4double worldXDimension = 1.*m;
106  G4double worldYDimension = 1.*m;
107  G4double worldZDimension = 1.*m;
108 
109  fWorld_solid = new G4Box( "WorldSolid",
110  worldXDimension,
111  worldYDimension,
112  worldZDimension );
113 
114  fWorld_logic = new G4LogicalVolume( fWorld_solid,
115  fAir,
116  "WorldLogical",
117  0, 0, 0 );
118 
119  fWorld_phys = new G4PVPlacement( 0,
120  G4ThreeVector(0,0,0),
121  "World",
122  fWorld_logic,
123  0,
124  false,
125  0 );
126 
127 #ifdef G4_DCMTK
128  ReadPhantomDataNew();
129  ConstructPhantomContainerNew();
130 #else
131  ReadPhantomData();
132  ConstructPhantomContainer();
133 #endif
134 
135  ConstructPhantom();
136  }
137  return fWorld_phys;
138 }
139 
140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
141 void DicomDetectorConstruction::InitialisationOfMaterials()
142 {
143  // Creating elements :
144  G4double z, a, density;
145  G4String name, symbol;
146 
147  G4Element* elC = new G4Element( name = "Carbon",
148  symbol = "C",
149  z = 6.0, a = 12.011 * g/mole );
150  G4Element* elH = new G4Element( name = "Hydrogen",
151  symbol = "H",
152  z = 1.0, a = 1.008 * g/mole );
153  G4Element* elN = new G4Element( name = "Nitrogen",
154  symbol = "N",
155  z = 7.0, a = 14.007 * g/mole );
156  G4Element* elO = new G4Element( name = "Oxygen",
157  symbol = "O",
158  z = 8.0, a = 16.00 * g/mole );
159  G4Element* elNa = new G4Element( name = "Sodium",
160  symbol = "Na",
161  z= 11.0, a = 22.98977* g/mole );
162  G4Element* elS = new G4Element( name = "Sulfur",
163  symbol = "S",
164  z = 16.0,a = 32.065* g/mole );
165  G4Element* elCl = new G4Element( name = "Chlorine",
166  symbol = "P",
167  z = 17.0, a = 35.453* g/mole );
168  G4Element* elK = new G4Element( name = "Potassium",
169  symbol = "P",
170  z = 19.0, a = 39.0983* g/mole );
171  G4Element* elP = new G4Element( name = "Phosphorus",
172  symbol = "P",
173  z = 15.0, a = 30.973976* g/mole );
174  G4Element* elFe = new G4Element( name = "Iron",
175  symbol = "Fe",
176  z = 26, a = 56.845* g/mole );
177  G4Element* elMg = new G4Element( name = "Magnesium",
178  symbol = "Mg",
179  z = 12.0, a = 24.3050* g/mole );
180  G4Element* elCa = new G4Element( name="Calcium",
181  symbol = "Ca",
182  z = 20.0, a = 40.078* g/mole );
183 
184  // Creating Materials :
185  G4int numberofElements;
186 
187  // Air
188  fAir = new G4Material( "Air",
189  1.290*mg/cm3,
190  numberofElements = 2 );
191  fAir->AddElement(elN, 0.7);
192  fAir->AddElement(elO, 0.3);
193 
194  // Lung Inhale
195  G4Material* lunginhale = new G4Material( "LungInhale",
196  density = 0.217*g/cm3,
197  numberofElements = 9);
198  lunginhale->AddElement(elH,0.103);
199  lunginhale->AddElement(elC,0.105);
200  lunginhale->AddElement(elN,0.031);
201  lunginhale->AddElement(elO,0.749);
202  lunginhale->AddElement(elNa,0.002);
203  lunginhale->AddElement(elP,0.002);
204  lunginhale->AddElement(elS,0.003);
205  lunginhale->AddElement(elCl,0.002);
206  lunginhale->AddElement(elK,0.003);
207 
208  // Lung exhale
209  G4Material* lungexhale = new G4Material( "LungExhale",
210  density = 0.508*g/cm3,
211  numberofElements = 9 );
212  lungexhale->AddElement(elH,0.103);
213  lungexhale->AddElement(elC,0.105);
214  lungexhale->AddElement(elN,0.031);
215  lungexhale->AddElement(elO,0.749);
216  lungexhale->AddElement(elNa,0.002);
217  lungexhale->AddElement(elP,0.002);
218  lungexhale->AddElement(elS,0.003);
219  lungexhale->AddElement(elCl,0.002);
220  lungexhale->AddElement(elK,0.003);
221 
222  // Adipose tissue
223  G4Material* adiposeTissue = new G4Material( "AdiposeTissue",
224  density = 0.967*g/cm3,
225  numberofElements = 7);
226  adiposeTissue->AddElement(elH,0.114);
227  adiposeTissue->AddElement(elC,0.598);
228  adiposeTissue->AddElement(elN,0.007);
229  adiposeTissue->AddElement(elO,0.278);
230  adiposeTissue->AddElement(elNa,0.001);
231  adiposeTissue->AddElement(elS,0.001);
232  adiposeTissue->AddElement(elCl,0.001);
233 
234  // Breast
235  G4Material* breast = new G4Material( "Breast",
236  density = 0.990*g/cm3,
237  numberofElements = 8 );
238  breast->AddElement(elH,0.109);
239  breast->AddElement(elC,0.506);
240  breast->AddElement(elN,0.023);
241  breast->AddElement(elO,0.358);
242  breast->AddElement(elNa,0.001);
243  breast->AddElement(elP,0.001);
244  breast->AddElement(elS,0.001);
245  breast->AddElement(elCl,0.001);
246 
247  // Water
248  G4Material* water = new G4Material( "Water",
249  density = 1.0*g/cm3,
250  numberofElements = 2 );
251  water->AddElement(elH,0.112);
252  water->AddElement(elO,0.888);
253 
254  // Muscle
255  G4Material* muscle = new G4Material( "Muscle",
256  density = 1.061*g/cm3,
257  numberofElements = 9 );
258  muscle->AddElement(elH,0.102);
259  muscle->AddElement(elC,0.143);
260  muscle->AddElement(elN,0.034);
261  muscle->AddElement(elO,0.710);
262  muscle->AddElement(elNa,0.001);
263  muscle->AddElement(elP,0.002);
264  muscle->AddElement(elS,0.003);
265  muscle->AddElement(elCl,0.001);
266  muscle->AddElement(elK,0.004);
267 
268  // Liver
269  G4Material* liver = new G4Material( "Liver",
270  density = 1.071*g/cm3,
271  numberofElements = 9);
272  liver->AddElement(elH,0.102);
273  liver->AddElement(elC,0.139);
274  liver->AddElement(elN,0.030);
275  liver->AddElement(elO,0.716);
276  liver->AddElement(elNa,0.002);
277  liver->AddElement(elP,0.003);
278  liver->AddElement(elS,0.003);
279  liver->AddElement(elCl,0.002);
280  liver->AddElement(elK,0.003);
281 
282  // Trabecular Bone
283  G4Material* trabecularBone = new G4Material( "TrabecularBone",
284  density = 1.159*g/cm3,
285  numberofElements = 12 );
286  trabecularBone->AddElement(elH,0.085);
287  trabecularBone->AddElement(elC,0.404);
288  trabecularBone->AddElement(elN,0.058);
289  trabecularBone->AddElement(elO,0.367);
290  trabecularBone->AddElement(elNa,0.001);
291  trabecularBone->AddElement(elMg,0.001);
292  trabecularBone->AddElement(elP,0.034);
293  trabecularBone->AddElement(elS,0.002);
294  trabecularBone->AddElement(elCl,0.002);
295  trabecularBone->AddElement(elK,0.001);
296  trabecularBone->AddElement(elCa,0.044);
297  trabecularBone->AddElement(elFe,0.001);
298 
299  // Dense Bone
300  G4Material* denseBone = new G4Material( "DenseBone",
301  density = 1.575*g/cm3,
302  numberofElements = 11 );
303  denseBone->AddElement(elH,0.056);
304  denseBone->AddElement(elC,0.235);
305  denseBone->AddElement(elN,0.050);
306  denseBone->AddElement(elO,0.434);
307  denseBone->AddElement(elNa,0.001);
308  denseBone->AddElement(elMg,0.001);
309  denseBone->AddElement(elP,0.072);
310  denseBone->AddElement(elS,0.003);
311  denseBone->AddElement(elCl,0.001);
312  denseBone->AddElement(elK,0.001);
313  denseBone->AddElement(elCa,0.146);
314 
315  //----- Put the materials in a vector
316  fOriginalMaterials.push_back(fAir); // rho = 0.00129
317  fOriginalMaterials.push_back(lunginhale); // rho = 0.217
318  fOriginalMaterials.push_back(lungexhale); // rho = 0.508
319  fOriginalMaterials.push_back(adiposeTissue); // rho = 0.967
320  fOriginalMaterials.push_back(breast ); // rho = 0.990
321  fOriginalMaterials.push_back(water); // rho = 1.018
322  fOriginalMaterials.push_back(muscle); // rho = 1.061
323  fOriginalMaterials.push_back(liver); // rho = 1.071
324  fOriginalMaterials.push_back(trabecularBone); // rho = 1.159
325  fOriginalMaterials.push_back(denseBone); // rho = 1.575
326 
327 }
328 
329 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
330 void DicomDetectorConstruction::ReadPhantomDataNew()
331 {
332 #ifdef G4_DCMTK
333  G4String fileName = DicomFileMgr::GetInstance()->GetFileOutName();
334 
335  std::ifstream fin(fileName);
336  std::vector<G4String> wl;
337  G4int nMaterials;
338  fin >> nMaterials;
339  G4String mateName;
340  G4int nmate;
341  for( G4int ii = 0; ii < nMaterials; ii++ ){
342  fin >> nmate;
343  fin >> mateName;
344  if( mateName[0] == '"' && mateName[mateName.length()-1] == '"' ) {
345  mateName = mateName.substr(1,mateName.length()-2);
346  }
347  G4cout << "GmReadPhantomG4Geometry::ReadPhantomData reading nmate " << ii << " = " << nmate
348  << " mate " << mateName << G4endl;
349  if( ii != nmate ) G4Exception("GmReadPhantomG4Geometry::ReadPhantomData",
350  "Wrong argument",
351  FatalErrorInArgument,
352  "Material number should be in increasing order: wrong material number ");
353 
354  G4Material* mate = 0;
355  const G4MaterialTable* matTab = G4Material::GetMaterialTable();
356  std::vector<G4Material*>::const_iterator matite;
357  for( matite = matTab->begin(); matite != matTab->end(); ++matite ) {
358  if( (*matite)->GetName() == mateName ) {
359  mate = *matite;
360  }
361  }
362  if( mate == 0 ) {
363  mate = G4NistManager::Instance()->FindOrBuildMaterial(mateName);
364  }
365  if( !mate ) G4Exception("GmReadPhantomG4Geometry::ReadPhantomData",
366  "Wrong argument",
367  FatalErrorInArgument,
368  ("Material not found" + mateName).c_str());
369  thePhantomMaterialsOriginal[nmate] = mate;
370  }
371 
372  fin >> fNVoxelX >> fNVoxelY >> fNVoxelZ;
373  G4cout << "GmReadPhantomG4Geometry::ReadPhantomData fNVoxel X/Y/Z " << fNVoxelX << " "
374  << fNVoxelY << " " << fNVoxelZ << G4endl;
375  fin >> fMinX >> fMaxX;
376  fin >> fMinY >> fMaxY;
377  fin >> fMinZ >> fMaxZ;
378  fVoxelHalfDimX = (fMaxX-fMinX)/fNVoxelX/2.;
379  fVoxelHalfDimY = (fMaxY-fMinY)/fNVoxelY/2.;
380  fVoxelHalfDimZ = (fMaxZ-fMinZ)/fNVoxelZ/2.;
381 #ifdef G4VERBOSE
382  G4cout << " Extension in X " << fMinX << " " << fMaxX << G4endl
383  << " Extension in Y " << fMinY << " " << fMaxY << G4endl
384  << " Extension in Z " << fMinZ << " " << fMaxZ << G4endl;
385 #endif
386 
387  fMateIDs = new size_t[fNVoxelX*fNVoxelY*fNVoxelZ];
388  for( G4int iz = 0; iz < fNVoxelZ; iz++ ) {
389  for( G4int iy = 0; iy < fNVoxelY; iy++ ) {
390  for( G4int ix = 0; ix < fNVoxelX; ix++ ) {
391  G4int mateID;
392  fin >> mateID;
393  G4int nnew = ix + (iy)*fNVoxelX + (iz)*fNVoxelX*fNVoxelY;
394  if( mateID < 0 || mateID >= nMaterials ) {
395  G4Exception("GmReadPhantomG4Geometry::ReadPhantomData",
396  "Wrong index in phantom file",
397  FatalException,
398  G4String("It should be between 0 and "
399  + G4UIcommand::ConvertToString(nMaterials-1)
400  + ", while it is "
401  + G4UIcommand::ConvertToString(mateID)).c_str());
402  }
403  fMateIDs[nnew] = mateID;
404  }
405  }
406  }
407 
408  ReadVoxelDensities( fin );
409 
410  fin.close();
411 #endif
412 
413 }
414 
415 
416 void DicomDetectorConstruction::ReadVoxelDensities( std::ifstream& fin )
417 {
418  G4String stemp;
419  std::map<G4int, std::pair<G4double,G4double> > densiMinMax;
420  std::map<G4int, std::pair<G4double,G4double> >::iterator mpite;
421  for( size_t ii = 0; ii < thePhantomMaterialsOriginal.size(); ii++ ){
422  densiMinMax[ii] = std::pair<G4double,G4double>(DBL_MAX,-DBL_MAX);
423  }
424 
425  char* part = getenv( "DICOM_CHANGE_MATERIAL_DENSITY" );
426  G4double densityDiff = -1.;
427  if( part ) densityDiff = G4UIcommand::ConvertToDouble(part);
428 
429  std::map<G4int,G4double> densityDiffs;
430  for( size_t ii = 0; ii < thePhantomMaterialsOriginal.size(); ii++ ){
431  densityDiffs[ii] = densityDiff; //currently all materials with same step
432  }
433  // densityDiffs[0] = 0.0001; //air
434 
435  //--- Calculate the average material density for each material/density bin
436  std::map< std::pair<G4Material*,G4int>, matInfo* > newMateDens;
437 
438  //---- Read the material densities
439  G4double dens;
440  for( G4int iz = 0; iz < fNVoxelZ; iz++ ) {
441  for( G4int iy = 0; iy < fNVoxelY; iy++ ) {
442  for( G4int ix = 0; ix < fNVoxelX; ix++ ) {
443  fin >> dens;
444  // G4cout << ix << " " << iy << " " << iz << " density " << dens << G4endl;
445  G4int copyNo = ix + (iy)*fNVoxelX + (iz)*fNVoxelX*fNVoxelY;
446 
447  if( densityDiff != -1. ) continue;
448 
449  //--- store the minimum and maximum density for each material (just for printing)
450  mpite = densiMinMax.find( fMateIDs[copyNo] );
451  if( dens < (*mpite).second.first ) (*mpite).second.first = dens;
452  if( dens > (*mpite).second.second ) (*mpite).second.second = dens;
453  //--- Get material from original list of material in file
454  int mateID = fMateIDs[copyNo];
455  std::map<G4int,G4Material*>::const_iterator imite =
456  thePhantomMaterialsOriginal.find(mateID);
457  // G4cout << copyNo << " mateID " << mateID << G4endl;
458  //--- Check if density is equal to the original material density
459  if( std::fabs(dens - (*imite).second->GetDensity()/CLHEP::g*CLHEP::cm3 ) < 1.e-9 ) continue;
460 
461  //--- Build material name with thePhantomMaterialsOriginal name + density
462  // float densityBin = densityDiffs[mateID] * (G4int(dens/densityDiffs[mateID])+0.5);
463  G4int densityBin = (G4int(dens/densityDiffs[mateID]));
464 
465  G4String mateName = (*imite).second->GetName()+G4UIcommand::ConvertToString(densityBin);
466  //--- Look if it is the first voxel with this material/densityBin
467  std::pair<G4Material*,G4int> matdens((*imite).second, densityBin );
468 
469  std::map< std::pair<G4Material*,G4int>, matInfo* >::iterator mppite =
470  newMateDens.find( matdens );
471  if( mppite != newMateDens.end() ){
472  matInfo* mi = (*mppite).second;
473  mi->fSumdens += dens;
474  mi->fNvoxels++;
475  fMateIDs[copyNo] = thePhantomMaterialsOriginal.size()-1 + mi->fId;
476  } else {
477  matInfo* mi = new matInfo;
478  mi->fSumdens = dens;
479  mi->fNvoxels = 1;
480  mi->fId = newMateDens.size()+1;
481  newMateDens[matdens] = mi;
482  fMateIDs[copyNo] = thePhantomMaterialsOriginal.size()-1 + mi->fId;
483  }
484  }
485  }
486  }
487 
488  if( densityDiff != -1. ) {
489  for( mpite = densiMinMax.begin(); mpite != densiMinMax.end(); mpite++ ){
490 #ifdef G4VERBOSE
491  G4cout << "DicomDetectorConstruction::ReadVoxelDensities ORIG MATERIALS DENSITY "
492  << (*mpite).first << " MIN " << (*mpite).second.first << " MAX "
493  << (*mpite).second.second << G4endl;
494 #endif
495  }
496  }
497 
498  //----- Build the list of phantom materials that go to Parameterisation
499  //--- Add original materials
500  std::map<G4int,G4Material*>::const_iterator mimite;
501  for( mimite = thePhantomMaterialsOriginal.begin(); mimite != thePhantomMaterialsOriginal.end();
502  mimite++ ){
503  fMaterials.push_back( (*mimite).second );
504  }
505  //
506  //---- Build and add new materials
507 std::map< std::pair<G4Material*,G4int>, matInfo* >::iterator mppite;
508  for( mppite= newMateDens.begin(); mppite != newMateDens.end(); mppite++ ){
509  G4double averdens = (*mppite).second->fSumdens/(*mppite).second->fNvoxels;
510  G4double saverdens = G4int(1000.001*averdens)/1000.;
511 #ifndef G4VERBOSE
512  G4cout << "DicomDetectorConstruction::ReadVoxelDensities AVER DENS " << averdens << " -> "
513  << saverdens << " -> " << G4int(1000*averdens) << " " << G4int(1000*averdens)/1000
514  << " " << G4int(1000*averdens)/1000. << G4endl;
515 #endif
516 
517  G4String mateName = ((*mppite).first).first->GetName() + "_"
518  + G4UIcommand::ConvertToString(saverdens);
519  fMaterials.push_back( BuildMaterialWithChangingDensity(
520  (*mppite).first.first, averdens, mateName ) );
521  }
522 
523 }
524 
525 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
526 void DicomDetectorConstruction::ReadPhantomData()
527 {
528 
529  G4String dataFile = "Data.dat";
530  std::ifstream finDF(dataFile.c_str());
531  G4String fname;
532  if(finDF.good() != 1 ) {
533  G4String descript = "Problem reading data file: "+dataFile;
534  G4Exception(" DicomDetectorConstruction::ReadPhantomData",
535  "",
536  FatalException,
537  descript);
538  }
539 
540  G4int compression;
541  finDF >> compression; // not used here
542 
543  finDF >> fNoFiles;
544  for(G4int i = 0; i < fNoFiles; i++ ) {
545  finDF >> fname;
546  //--- Read one data file
547  fname += ".g4dcm";
548  ReadPhantomDataFile(fname);
549  }
550 
551  //----- Merge data headers
552  MergeZSliceHeaders();
553 
554  finDF.close();
555 
556 }
557 
558 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
559 void DicomDetectorConstruction::ReadPhantomDataFile(const G4String& fname)
560 {
561  G4cout << " DicomDetectorConstruction::ReadPhantomDataFile opening file "
562  << fname << G4endl; //GDEB
563 #ifdef G4VERBOSE
564  G4cout << " DicomDetectorConstruction::ReadPhantomDataFile opening file "
565  << fname << G4endl;
566 #endif
567  std::ifstream fin(fname.c_str(), std::ios_base::in);
568  if( !fin.is_open() ) {
569  G4Exception("DicomDetectorConstruction::ReadPhantomDataFile",
570  "",
571  FatalErrorInArgument,
572  G4String("File not found " + fname ).c_str());
573  }
574  //----- Define density differences (maximum density difference to create
575  // a new material)
576  char* part = getenv( "DICOM_CHANGE_MATERIAL_DENSITY" );
577  G4double densityDiff = -1.;
578  if( part ) densityDiff = G4UIcommand::ConvertToDouble(part);
579  if( densityDiff != -1. ) {
580  for( unsigned int ii = 0; ii < fOriginalMaterials.size(); ii++ ){
581  fDensityDiffs[ii] = densityDiff; //currently all materials with
582  // same difference
583  }
584  }else {
585  if( fMaterials.size() == 0 ) { // do it only for first slice
586  for( unsigned int ii = 0; ii < fOriginalMaterials.size(); ii++ ){
587  fMaterials.push_back( fOriginalMaterials[ii] );
588  }
589  }
590  }
591 
592  //----- Read data header
593  DicomPhantomZSliceHeader* sliceHeader = new DicomPhantomZSliceHeader( fin );
594  fZSliceHeaders.push_back( sliceHeader );
595 
596  //----- Read material indices
597  G4int nVoxels = sliceHeader->GetNoVoxels();
598 
599  //--- If first slice, initiliaze fMateIDs
600  if( fZSliceHeaders.size() == 1 ) {
601  //fMateIDs = new unsigned int[fNoFiles*nVoxels];
602  fMateIDs = new size_t[fNoFiles*nVoxels];
603 
604  }
605 
606  unsigned int mateID;
607  // number of voxels from previously read slices
608  G4int voxelCopyNo = (fZSliceHeaders.size()-1)*nVoxels;
609  for( G4int ii = 0; ii < nVoxels; ii++, voxelCopyNo++ ){
610  fin >> mateID;
611  fMateIDs[voxelCopyNo] = mateID;
612  }
613 
614  //----- Read material densities and build new materials if two voxels have
615  // same material but its density is in a different density interval
616  // (size of density intervals defined by densityDiff)
617  G4double density;
618  // number of voxels from previously read slices
619  voxelCopyNo = (fZSliceHeaders.size()-1)*nVoxels;
620  for( G4int ii = 0; ii < nVoxels; ii++, voxelCopyNo++ ){
621  fin >> density;
622 
623  //-- Get material from list of original materials
624  mateID = fMateIDs[voxelCopyNo];
625  G4Material* mateOrig = fOriginalMaterials[mateID];
626 
627  //-- Get density bin: middle point of the bin in which the current
628  // density is included
629  G4String newMateName = mateOrig->GetName();
630  float densityBin = 0.;
631  if( densityDiff != -1.) {
632  densityBin = fDensityDiffs[mateID] *
633  (G4int(density/fDensityDiffs[mateID])+0.5);
634  //-- Build the new material name
635  newMateName += G4UIcommand::ConvertToString(densityBin);
636  }
637 
638  //-- Look if a material with this name is already created
639  // (because a previous voxel was already in this density bin)
640  unsigned int im;
641  for( im = 0; im < fMaterials.size(); im++ ){
642  if( fMaterials[im]->GetName() == newMateName ) {
643  break;
644  }
645  }
646  //-- If material is already created use index of this material
647  if( im != fMaterials.size() ) {
648  fMateIDs[voxelCopyNo] = im;
649  //-- else, create the material
650  } else {
651  if( densityDiff != -1.) {
652  fMaterials.push_back( BuildMaterialWithChangingDensity( mateOrig,
653  densityBin, newMateName ) );
654  fMateIDs[voxelCopyNo] = fMaterials.size()-1;
655  } else {
656  G4cerr << " im " << im << " < " << fMaterials.size() << " name "
657  << newMateName << G4endl;
658  G4Exception("DicomDetectorConstruction::ReadPhantomDataFile",
659  "",
660  FatalErrorInArgument,
661  "Wrong index in material"); //it should never reach here
662  }
663  }
664  }
665 
666 }
667 
668 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
669 void DicomDetectorConstruction::MergeZSliceHeaders()
670 {
671  //----- Images must have the same dimension ...
672  fZSliceHeaderMerged = new DicomPhantomZSliceHeader( *fZSliceHeaders[0] );
673  for( unsigned int ii = 1; ii < fZSliceHeaders.size(); ii++ ) {
674  *fZSliceHeaderMerged += *fZSliceHeaders[ii];
675  };
676 
677 }
678 
679 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
680 G4Material* DicomDetectorConstruction::BuildMaterialWithChangingDensity(
681  const G4Material* origMate, float density, G4String newMateName )
682 {
683  //----- Copy original material, but with new density
684  G4int nelem = origMate->GetNumberOfElements();
685  G4Material* mate = new G4Material( newMateName, density*g/cm3, nelem,
686  kStateUndefined, STP_Temperature );
687 
688  for( G4int ii = 0; ii < nelem; ii++ ){
689  G4double frac = origMate->GetFractionVector()[ii];
690  G4Element* elem = const_cast<G4Element*>(origMate->GetElement(ii));
691  mate->AddElement( elem, frac );
692  }
693 
694  return mate;
695 }
696 
697 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
698 void DicomDetectorConstruction::ConstructPhantomContainer()
699 {
700  //---- Extract number of voxels and voxel dimensions
701  fNVoxelX = fZSliceHeaderMerged->GetNoVoxelX();
702  fNVoxelY = fZSliceHeaderMerged->GetNoVoxelY();
703  fNVoxelZ = fZSliceHeaderMerged->GetNoVoxelZ();
704 
705  fVoxelHalfDimX = fZSliceHeaderMerged->GetVoxelHalfX();
706  fVoxelHalfDimY = fZSliceHeaderMerged->GetVoxelHalfY();
707  fVoxelHalfDimZ = fZSliceHeaderMerged->GetVoxelHalfZ();
708 #ifdef G4VERBOSE
709  G4cout << " fNVoxelX " << fNVoxelX << " fVoxelHalfDimX " << fVoxelHalfDimX
710  <<G4endl;
711  G4cout << " fNVoxelY " << fNVoxelY << " fVoxelHalfDimY " << fVoxelHalfDimY
712  <<G4endl;
713  G4cout << " fNVoxelZ " << fNVoxelZ << " fVoxelHalfDimZ " << fVoxelHalfDimZ
714  <<G4endl;
715  G4cout << " totalPixels " << fNVoxelX*fNVoxelY*fNVoxelZ << G4endl;
716 #endif
717 
718  //----- Define the volume that contains all the voxels
719  fContainer_solid = new G4Box("phantomContainer",fNVoxelX*fVoxelHalfDimX,
720  fNVoxelY*fVoxelHalfDimY,
721  fNVoxelZ*fVoxelHalfDimZ);
722  fContainer_logic =
723  new G4LogicalVolume( fContainer_solid,
724  //the material is not important, it will be fully filled by the voxels
725  fMaterials[0],
726  "phantomContainer",
727  0, 0, 0 );
728  //--- Place it on the world
729  G4double fOffsetX = (fZSliceHeaderMerged->GetMaxX() +
730  fZSliceHeaderMerged->GetMinX() ) /2.;
731  G4double fOffsetY = (fZSliceHeaderMerged->GetMaxY() +
732  fZSliceHeaderMerged->GetMinY() ) /2.;
733  G4double fOffsetZ = (fZSliceHeaderMerged->GetMaxZ() +
734  fZSliceHeaderMerged->GetMinZ() ) /2.;
735  G4ThreeVector posCentreVoxels(fOffsetX,fOffsetY,fOffsetZ);
736 #ifdef G4VERBOSE
737  G4cout << " placing voxel container volume at " << posCentreVoxels << G4endl;
738 #endif
739  fContainer_phys =
740  new G4PVPlacement(0, // rotation
741  posCentreVoxels,
742  fContainer_logic, // The logic volume
743  "phantomContainer", // Name
744  fWorld_logic, // Mother
745  false, // No op. bool.
746  1); // Copy number
747 
748  //fContainer_logic->SetVisAttributes(new G4VisAttributes(G4Colour(1.,0.,0.)));
749 }
750 
751 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
752 void DicomDetectorConstruction::ConstructPhantomContainerNew()
753 {
754 #ifdef G4_DCMTK
755  //---- Extract number of voxels and voxel dimensions
756 #ifdef G4VERBOSE
757  G4cout << " fNVoxelX " << fNVoxelX << " fVoxelHalfDimX " << fVoxelHalfDimX
758  <<G4endl;
759  G4cout << " fNVoxelY " << fNVoxelY << " fVoxelHalfDimY " << fVoxelHalfDimY
760  <<G4endl;
761  G4cout << " fNVoxelZ " << fNVoxelZ << " fVoxelHalfDimZ " << fVoxelHalfDimZ
762  <<G4endl;
763  G4cout << " totalPixels " << fNVoxelX*fNVoxelY*fNVoxelZ << G4endl;
764 #endif
765 
766  //----- Define the volume that contains all the voxels
767  fContainer_solid = new G4Box("phantomContainer",fNVoxelX*fVoxelHalfDimX,
768  fNVoxelY*fVoxelHalfDimY,
769  fNVoxelZ*fVoxelHalfDimZ);
770  fContainer_logic =
771  new G4LogicalVolume( fContainer_solid,
772  //the material is not important, it will be fully filled by the voxels
773  fMaterials[0],
774  "phantomContainer",
775  0, 0, 0 );
776 
777  G4ThreeVector posCentreVoxels((fMinX+fMaxX)/2.,(fMinY+fMaxY)/2.,(fMinZ+fMaxZ)/2.);
778 #ifdef G4VERBOSE
779  G4cout << " placing voxel container volume at " << posCentreVoxels << G4endl;
780 #endif
781  fContainer_phys =
782  new G4PVPlacement(0, // rotation
783  posCentreVoxels,
784  fContainer_logic, // The logic volume
785  "phantomContainer", // Name
786  fWorld_logic, // Mother
787  false, // No op. bool.
788  1); // Copy number
789 
790  //fContainer_logic->SetVisAttributes(new G4VisAttributes(G4Colour(1.,0.,0.)));
791 #endif
792 }
793 
794 #include "G4SDManager.hh"
795 #include "G4MultiFunctionalDetector.hh"
796 #include "G4PSDoseDeposit.hh"
797 #include "G4PSDoseDeposit3D.hh"
798 
799 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
800 void DicomDetectorConstruction::SetScorer(G4LogicalVolume* voxel_logic)
801 {
802 
803 #ifdef G4VERBOSE
804  G4cout << "\t SET SCORER : " << voxel_logic->GetName() << G4endl;
805 #endif
806 
807  fScorers.insert(voxel_logic);
808 
809 }
810 
811 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
812 
813 void DicomDetectorConstruction::ConstructSDandField()
814 {
815 
816 #ifdef G4VERBOSE
817  G4cout << "\t CONSTRUCT SD AND FIELD" << G4endl;
818 #endif
819 
820  //G4SDManager* SDman = G4SDManager::GetSDMpointer();
821 
822  //SDman->SetVerboseLevel(1);
823 
824  //
825  // Sensitive Detector Name
826  G4String concreteSDname = "phantomSD";
827  std::vector<G4String> scorer_names;
828  scorer_names.push_back(concreteSDname);
829  //------------------------
830  // MultiFunctionalDetector
831  //------------------------
832  //
833  // Define MultiFunctionalDetector with name.
834  // declare MFDet as a MultiFunctionalDetector scorer
835  G4MultiFunctionalDetector* MFDet =
836  new G4MultiFunctionalDetector(concreteSDname);
837  G4SDManager::GetSDMpointer()->AddNewDetector( MFDet );
838  //G4VPrimitiveScorer* dosedep = new G4PSDoseDeposit("DoseDeposit");
839  G4VPrimitiveScorer* dosedep =
840  new G4PSDoseDeposit3D("DoseDeposit", fNVoxelX, fNVoxelY, fNVoxelZ);
841  MFDet->RegisterPrimitive(dosedep);
842 
843  for(std::set<G4LogicalVolume*>::iterator ite = fScorers.begin();
844  ite != fScorers.end(); ++ite) {
845  SetSensitiveDetector(*ite, MFDet);
846  }
847 
848  /*if(DicomRunAction::Instance()->GetDicomRun()) {
849  DicomRunAction::Instance()->GetDicomRun()->ConstructMFD(scorer_names);
850  }*/
851 
852 }
CLHEP::m
static constexpr double m
Definition: SystemOfUnits.h:109
name
const XML_Char * name
Definition: expat.h:151
G4MultiFunctionalDetector::RegisterPrimitive
G4bool RegisterPrimitive(G4VPrimitiveScorer *)
Definition: G4MultiFunctionalDetector.cc:57
DicomDetectorConstruction::SetScorer
void SetScorer(G4LogicalVolume *voxel_logic)
Definition: DicomDetectorConstruction.cc:800
G4NistManager::FindOrBuildMaterial
G4Material * FindOrBuildMaterial(const G4String &name, G4bool isotopes=true, G4bool warning=false)
Definition: G4NistManager.hh:474
DicomDetectorConstruction.hh
Definition of the DicomDetectorConstruction class.
DicomDetectorConstruction::fScorers
std::set< G4LogicalVolume * > fScorers
Definition: DicomDetectorConstruction.hh:144
G4MultiFunctionalDetector.hh
kStateUndefined
Definition: G4Material.hh:114
mg
static constexpr double mg
Definition: G4SIunits.hh:184
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4Element.hh
a
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
G4Box
Definition: G4Box.hh:64
DicomRunAction.hh
Definition of the DicomRunAction class.
G4NistManager.hh
G4Material::GetName
const G4String & GetName() const
Definition: G4Material.hh:178
DicomDetectorConstruction::fZSliceHeaders
std::vector< DicomPhantomZSliceHeader * > fZSliceHeaders
Definition: DicomDetectorConstruction.hh:129
G4Material::GetMaterialTable
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:587
G4MaterialTable
std::vector< G4Material * > G4MaterialTable
Definition: G4MaterialTable.hh:41
first
G4int first
Definition: G4ParticleHPJENDLHEData.cc:280
G4UIcommand::ConvertToString
static G4String ConvertToString(G4bool boolVal)
Definition: G4UIcommand.cc:372
G4Box.hh
DicomDetectorConstruction::Construct
virtual G4VPhysicalVolume * Construct()
Definition: DicomDetectorConstruction.cc:98
G4PSDoseDeposit.hh
python.hepunit.STP_Temperature
float STP_Temperature
Definition: hepunit.py:302
CLHEP::mg
static constexpr double mg
Definition: SystemOfUnits.h:191
G4Material::GetElement
const G4Element * GetElement(G4int iel) const
Definition: G4Material.hh:202
G4int
int G4int
Definition: G4Types.hh:78
G4NistManager::Instance
static G4NistManager * Instance()
Definition: G4NistManager.cc:68
DicomRun.hh
Definition of the DicomRun class.
CLHEP::g
static constexpr double g
Definition: SystemOfUnits.h:190
DicomDetectorConstruction::fZSliceHeaderMerged
DicomPhantomZSliceHeader * fZSliceHeaderMerged
Definition: DicomDetectorConstruction.hh:131
DicomDetectorConstruction::ReadVoxelDensities
void ReadVoxelDensities(std::ifstream &fin)
Definition: DicomDetectorConstruction.cc:416
g
function g(Y1, Y2, PT2)
Definition: hijing1.383.f:5205
G4VisAttributes.hh
DicomFileMgr.hh
DicomDetectorConstruction::BuildMaterialWithChangingDensity
G4Material * BuildMaterialWithChangingDensity(const G4Material *origMate, float density, G4String newMateName)
Definition: DicomDetectorConstruction.cc:680
G4cout
G4GLOB_DLL std::ostream G4cout
m
static constexpr double m
Definition: G4SIunits.hh:129
G4UIcommand.hh
G4SDManager.hh
G4Material.hh
G4UIcommand::ConvertToDouble
static G4double ConvertToDouble(const char *st)
Definition: G4UIcommand.cc:455
DicomFileMgr::GetFileOutName
G4String GetFileOutName() const
Definition: DicomFileMgr.hh:94
cm3
static constexpr double cm3
Definition: G4SIunits.hh:121
G4VUserDetectorConstruction::SetSensitiveDetector
void SetSensitiveDetector(const G4String &logVolName, G4VSensitiveDetector *aSD, G4bool multi=false)
Definition: G4VUserDetectorConstruction.cc:206
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4PhysicalConstants.hh
G4SDManager::AddNewDetector
void AddNewDetector(G4VSensitiveDetector *aSD)
Definition: G4SDManager.cc:71
G4LogicalVolume.hh
FatalErrorInArgument
Definition: G4ExceptionSeverity.hh:61
G4SDManager::GetSDMpointer
static G4SDManager * GetSDMpointer()
Definition: G4SDManager.cc:40
FatalException
Definition: G4ExceptionSeverity.hh:60
test.z
tuple z
Definition: test.py:28
G4endl
#define G4endl
Definition: G4ios.hh:61
G4Material::AddElement
void AddElement(G4Element *element, G4int nAtoms)
Definition: G4Material.cc:362
G4Material::GetNumberOfElements
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
test.fname
string fname
Definition: test.py:308
CLHEP::mole
static constexpr double mole
Definition: SystemOfUnits.h:266
DicomPhantomZSliceHeader.hh
Definition of the DicomPhantomZSliceHeader class.
G4VPhysicalVolume.hh
G4LogicalVolume::GetName
const G4String & GetName() const
G4double
double G4double
Definition: G4Types.hh:76
DicomDetectorConstruction::fMaterials
std::vector< G4Material * > fMaterials
Definition: DicomDetectorConstruction.hh:120
DicomDetectorConstruction::fOriginalMaterials
std::vector< G4Material * > fOriginalMaterials
Definition: DicomDetectorConstruction.hh:119
DicomDetectorConstruction::ReadPhantomDataFile
void ReadPhantomDataFile(const G4String &fname)
Definition: DicomDetectorConstruction.cc:559
G4PSDoseDeposit3D.hh
G4PVPlacement.hh
G4Material::GetFractionVector
const G4double * GetFractionVector() const
Definition: G4Material.hh:194
SystemOfUnits.h
DBL_MAX
#define DBL_MAX
Definition: templates.hh:83
mole
static constexpr double mole
Definition: G4SIunits.hh:286
DicomDetectorConstruction::ConstructPhantom
virtual void ConstructPhantom()=0
DicomFileMgr::GetInstance
static DicomFileMgr * GetInstance()
Definition: DicomFileMgr.cc:41
G4cerr
G4GLOB_DLL std::ostream G4cerr
DicomDetectorConstruction::fDensityDiffs
std::map< G4int, G4double > fDensityDiffs
Definition: DicomDetectorConstruction.hh:126
DicomDetectorConstruction::thePhantomMaterialsOriginal
std::map< G4int, G4Material * > thePhantomMaterialsOriginal
Definition: DicomDetectorConstruction.hh:139
CLHEP::cm3
static constexpr double cm3
Definition: SystemOfUnits.h:101