Geant4_10
DicomPartialDetectorConstruction.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: DicomPartialDetectorConstruction.cc 73076 2013-08-16 07:45:30Z gcosmo $
27 //
30 //
31 
32 #include "globals.hh"
33 
35 
36 #include "G4Box.hh"
37 #include "G4LogicalVolume.hh"
38 #include "G4VPhysicalVolume.hh"
39 #include "G4PVPlacement.hh"
40 #include "G4PVParameterised.hh"
41 #include "G4Material.hh"
42 #include "G4Element.hh"
43 #include "G4VisAttributes.hh"
44 #include "G4Colour.hh"
45 #include "G4ios.hh"
47 #include "G4NistManager.hh"
48 #include "G4UIcommand.hh"
49 #include "G4Tubs.hh"
50 
51 #include "G4tgbVolumeMgr.hh"
52 #include "G4tgbMaterialMgr.hh"
53 #include "G4tgrMessenger.hh"
54 #include "G4tgrMessenger.hh"
55 
56 #include "G4SystemOfUnits.hh"
57 
58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
61  fPartialPhantomParam(0),
62  fNVoxels(0),
63  fDimX(0),
64  fDimY(0),
65  fDimZ(0),
66  fOffsetX(0),
67  fOffsetY(0),
68  fOffsetZ(0),
69  fMateIDs(0)
70 {
71 }
72 
73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
75 {
76 }
77 
78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
80 {
82 
83  //----- Build world
84  G4double worldXDimension = 1.*m;
85  G4double worldYDimension = 1.*m;
86  G4double worldZDimension = 1.*m;
87 
88  G4Box* world_box = new G4Box( "WorldSolid",
89  worldXDimension,
90  worldYDimension,
91  worldZDimension );
92 
93  fWorld_logic = new G4LogicalVolume( world_box,
94  fAir,
95  "WorldLogical",
96  0, 0, 0 );
97 
98  G4VPhysicalVolume* world_phys = new G4PVPlacement( 0,
99  G4ThreeVector(0,0,0),
100  "World",
101  fWorld_logic,
102  0,
103  false,
104  0 );
105 
106  ReadPhantomData();
107 
108  ConstructPhantomContainer();
109  ConstructPhantom();
110 
111  return world_phys;
112 }
113 
114 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
115 void DicomPartialDetectorConstruction::ReadPhantomData()
116 {
117  G4String dataFile = "Data.dat";
118  std::ifstream finDF(dataFile.c_str());
119  G4String fname;
120  if(finDF.good() != 1 ) {
121  G4Exception(" DicomPartialDetectorConstruction::ReadPhantomData",
122  "",
124  " Problem reading data file: Data.dat");
125  }
126 
127  G4int compression;
128  finDF >> compression; // not used here
129 
130  G4int numFiles;
131  finDF >> numFiles; // only 1 file supported for the moment
132  if( numFiles != 1 ) {
133  G4Exception("DicomPartialDetectorConstruction::ReadPhantomData",
134  "",
136  "More than 1 DICOM file is not supported");
137  }
138  for(G4int i = 0; i < numFiles; i++ ) {
139  finDF >> fname;
140  //--- Read one data file
141  fname += ".g4pdcm";
142  ReadPhantomDataFile(fname);
143  }
144 
145  finDF.close();
146 
147 }
148 
149 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
150 void DicomPartialDetectorConstruction::ReadPhantomDataFile(const G4String& fname)
151 {
152  // G4String filename = "phantom.g4pdcm";
153 #ifdef G4VERBOSE
154  G4cout << " DicomDetectorConstruction::ReadPhantomDataFile opening file " << fname << G4endl;
155 #endif
156  std::ifstream fin(fname.c_str(), std::ios_base::in);
157  if( !fin.is_open() ) {
158  G4Exception("DicomPartialDetectorConstruction::ReadPhantomDataFile",
159  "",
161  G4String("File not found " + fname).c_str());
162  }
163  G4int nMaterials;
164  fin >> nMaterials;
165  G4String stemp;
166  G4int nmate;
167  for( G4int ii = 0; ii < nMaterials; ii++ ){
168  fin >> nmate >> stemp;
169  G4cout << "DicomPartialDetectorConstruction::ReadPhantomData reading nmate "
170  << ii << " = " << nmate << " mate " << stemp << G4endl;
171  if( ii != nmate ) G4Exception("DicomPartialDetectorConstruction::ReadPhantomData",
172  "",
174  "Material number should be in increasing order: \
175  wrong material number ");
176 
177  }
178 
179  fin >> fNVoxelX >> fNVoxelY >> fNVoxelZ;
180  G4cout << "DicomPartialDetectorConstruction::ReadPhantomData nVoxel X/Y/Z "
181  << fNVoxelX << " " << fNVoxelY << " " << fNVoxelZ << G4endl;
182  fin >> fOffsetX >> fDimX;
183  fDimX = (fDimX - fOffsetX)/fNVoxelX;
184  fin >> fOffsetY >> fDimY;
185  fDimY = (fDimY - fOffsetY)/fNVoxelY;
186  fin >> fOffsetZ >> fDimZ;
187  fDimZ = (fDimZ - fOffsetZ)/fNVoxelZ;
188  G4cout << "DicomPartialDetectorConstruction::ReadPhantomData voxelDimX "
189  << fDimX << " fOffsetX " << fOffsetX << G4endl;
190  G4cout << "DicomPartialDetectorConstruction::ReadPhantomData voxelDimY "
191  << fDimY << " fOffsetY " << fOffsetY << G4endl;
192  G4cout << "DicomPartialDetectorConstruction::ReadPhantomData voxelDimZ "
193  << fDimZ << " fOffsetZ " << fOffsetZ << G4endl;
194 
195  //--- Read voxels that are filled
196  fNVoxels = 0;
197  // G4bool* isFilled = new G4bool[fNVoxelX*fNVoxelY*fNVoxelZ];
198  // fFilledIDs = new size_t[fNVoxelZ*fNVoxelY+1];
199  //? fFilledIDs.insert(0);
200  int ifxmin1, ifxmax1;
201  for( G4int iz = 0; iz < fNVoxelZ; iz++ ) {
202  std::map< G4int, G4int > ifmin, ifmax;
203  for( G4int iy = 0; iy < fNVoxelY; iy++ ) {
204  fin >> ifxmin1 >> ifxmax1;
205  // check coherence ...
206 
207  ifmin[iy] = ifxmin1;
208  ifmax[iy] = ifxmax1;
209  G4int ifxdiff = ifxmax1-ifxmin1+1;
210  if( ifxmax1 == -1 && ifxmin1 == -1 ) ifxdiff = 0;
211  // fFilledIDs.insert(std::multimap<G4int,size_t>::value_type(ifxdiff+fNVoxels, ifxmin1));
212  fFilledIDs.insert(std::pair<size_t,size_t>(ifxdiff+fNVoxels-1, ifxmin1));
213  G4cout << "DicomPartialDetectorConstruction::ReadPhantomData insert FilledIDs "
214  << ifxdiff+fNVoxels-1 << " min " << ifxmin1 << " N= " << fFilledIDs.size() << G4endl;
215  //filledIDs[iz*fNVoxelY+iy+1] = ifxmax1-ifxmin1+1;
216  for( G4int ix = 0; ix < fNVoxelX; ix++ ) {
217  // G4cout << "REGNAV DicomPartialDetectorConstruction::ReadPhantomData checking
218  // to add voxel iz " << iz << " iy " << iy << " ix " << ix << " min= "
219  // << ifxmin1 << " max= " << ifxmax1 << G4endl;
220  // G4int copyNo = ix + (iy)*fNVoxelX + (iz)*fNVoxelX*fNVoxelY;
221  if( ix >= G4int(ifxmin1) && ix <= G4int(ifxmax1) ) {
222  // isFilled[copyNo] = true;
223  fNVoxels++;
224  } else {
225  // isFilled[copyNo] = false;
226  }
227  }
228  }
229  fFilledMins[iz] = ifmin;
230  fFilledMaxs[iz] = ifmax;
231  }
232 
233  //--- Read material IDs
234  G4int mateID1;
235  fMateIDs = new size_t[fNVoxelX*fNVoxelY*fNVoxelZ];
236  G4int copyNo = 0;
237  for( G4int iz = 0; iz < fNVoxelZ; iz++ ) {
238  std::map< G4int, G4int > ifmin = fFilledMins[iz];
239  std::map< G4int, G4int > ifmax = fFilledMaxs[iz];
240  for( G4int iy = 0; iy < fNVoxelY; iy++ ) {
241  ifxmin1 = ifmin[iy];
242  ifxmax1 = ifmax[iy];
243  for( G4int ix = 0; ix < fNVoxelX; ix++ ) {
244  if( ix >= G4int(ifxmin1) && ix <= G4int(ifxmax1) ) {
245  fin >> mateID1;
246  if( mateID1 < 0 || mateID1 >= nMaterials ) {
247  G4Exception("DicomPartialDetectorConstruction::ReadPhantomData",
248  "",
250  G4String("Wrong index in phantom file: It should be between 0 and "
251  + G4UIcommand::ConvertToString(nMaterials-1)
252  + ", while it is "
253  + G4UIcommand::ConvertToString(mateID1)).c_str());
254  }
255  fMateIDs[copyNo] = mateID1;
256  copyNo++;
257  }
258  }
259  }
260  }
261 
262  ReadVoxelDensitiesPartial( fin );
263 
264  fin.close();
265 }
266 
267 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
268 void DicomPartialDetectorConstruction::ReadVoxelDensitiesPartial( std::ifstream& fin )
269 {
270  std::map<G4int, std::pair<G4double,G4double> > densiMinMax;
271  std::map<G4int, std::pair<G4double,G4double> >::iterator mpite;
272  for( size_t ii = 0; ii < fOriginalMaterials.size(); ii++ ){
273  densiMinMax[ii] = std::pair<G4double,G4double>(DBL_MAX,-DBL_MAX);
274  }
275 
276  //----- Define density differences (maximum density difference to create a new material)
277  char* part = getenv( "DICOM_CHANGE_MATERIAL_DENSITY" );
278  G4double densityDiff = -1.;
279  if( part ) densityDiff = G4UIcommand::ConvertToDouble(part);
280  std::map<G4int,G4double> densityDiffs;
281  if( densityDiff != -1. ) {
282  for( size_t ii = 0; ii < fOriginalMaterials.size(); ii++ ){
283  densityDiffs[ii] = densityDiff; //currently all materials with same step
284  }
285  }else {
286  if( fMaterials.size() == 0 ) { // do it only for first slice
287  for( unsigned int ii = 0; ii < fOriginalMaterials.size(); ii++ ){
288  fMaterials.push_back( fOriginalMaterials[ii] );
289  }
290  }
291  }
292 
293  // densityDiffs[0] = 0.0001; //fAir
294 
295  //--- Calculate the average material density for each material/density bin
296  std::map< std::pair<G4Material*,G4int>, matInfo* > newMateDens;
297 
298  G4double dens1;
299  size_t ifxmin1, ifxmax1;
300 
301  //---- Read the material densities
302  G4int copyNo = 0;
303  for( G4int iz = 0; iz < fNVoxelZ; iz++ ) {
304  std::map< G4int, G4int > ifmin = fFilledMins[iz];
305  std::map< G4int, G4int > ifmax = fFilledMaxs[iz];
306  for( G4int iy = 0; iy < fNVoxelY; iy++ ) {
307  ifxmin1 = ifmin[iy];
308  ifxmax1 = ifmax[iy];
309  for( G4int ix = 0; ix < fNVoxelX; ix++ ) {
310  if( ix >= G4int(ifxmin1) && ix <= G4int(ifxmax1) ) {
311  fin >> dens1;
312  // G4cout << ix << " " << iy << " " << iz << " filling fMateIDs "
313  // << copyNo << " = " << atoi(stemp.c_str())-1 << " " << stemp << G4endl;
314 
315  //--- store the minimum and maximum density for each material (just for printing)
316  mpite = densiMinMax.find( fMateIDs[copyNo] );
317  if( dens1 < (*mpite).second.first ) (*mpite).second.first = dens1;
318  if( dens1 > (*mpite).second.second ) (*mpite).second.second = dens1;
319 
320  //--- Get material from original list of material in file
321  int mateID = fMateIDs[copyNo];
322  G4Material* mate = fOriginalMaterials[mateID];
323  // G4cout << copyNo << " mateID " << mateID << G4endl;
324  //--- Check if density is equal to the original material density
325  if( std::fabs(dens1 - mate->GetDensity()/g*cm3 ) < 1.e-9 ) continue;
326 
327  //--- Build material name with fOriginalMaterials name + density
328  // float densityBin = densityDiffs[mateID] *
329  // (G4int(dens1/densityDiffs[mateID])+0.5);
330  G4String newMateName = mate->GetName();
331 
332  G4int densityBin = 0;
333  // G4cout << " densityBin " << densityBin << " "
334  // << dens1 << " " <<densityDiffs[mateID] << G4endl;
335 
336  if( densityDiff != -1.) {
337  densityBin = (G4int(dens1/densityDiffs[mateID]));
338  newMateName = mate->GetName()+G4UIcommand::ConvertToString(densityBin);
339  }
340 
341  //--- Look if it is the first voxel with this material/densityBin
342  std::pair<G4Material*,G4int> matdens(mate, densityBin );
343 
344  std::map< std::pair<G4Material*,G4int>, matInfo* >::iterator mppite
345  = newMateDens.find( matdens );
346  if( mppite != newMateDens.end() ){
347  matInfo* mi = (*mppite).second;
348  mi->sumdens += dens1;
349  mi->nvoxels++;
350  fMateIDs[copyNo] = fOriginalMaterials.size()-1 + mi->id;
351  // G4cout << copyNo << " mat new again "
352  // << fOriginalMaterials.size()-1 + mi->id << " " << mi->id << G4endl;
353  } else {
354  matInfo* mi = new matInfo;
355  mi->sumdens = dens1;
356  mi->nvoxels = 1;
357  mi->id = newMateDens.size()+1;
358  newMateDens[matdens] = mi;
359  fMateIDs[copyNo] = fOriginalMaterials.size()-1 + mi->id;
360  // G4cout << copyNo << " mat new first "
361  // << fOriginalMaterials.size()-1 + mi->id << G4endl;
362  }
363  copyNo++;
364  // G4cout << ix << " " << iy << " " << iz
365  // << " filling fMateIDs " << copyNo << " = " << atoi(cid)-1 << G4endl;
366  // fMateIDs[copyNo] = atoi(cid)-1;
367  }
368  }
369  }
370  }
371 
372  //----- Build the list of phantom materials that go to Parameterisation
373  //--- Add original materials
374  std::vector<G4Material*>::const_iterator mimite;
375  for( mimite = fOriginalMaterials.begin(); mimite != fOriginalMaterials.end(); mimite++ ){
376  fPhantomMaterials.push_back( (*mimite) );
377  }
378  //
379  //---- Build and add new materials
380  std::map< std::pair<G4Material*,G4int>, matInfo* >::iterator mppite;
381  for( mppite= newMateDens.begin(); mppite != newMateDens.end(); mppite++ ){
382  G4double averdens = (*mppite).second->sumdens/(*mppite).second->nvoxels;
383  G4double saverdens = G4int(1000.001*averdens)/1000.;
384  G4cout << "GmReadPhantomGeometry::ReadVoxelDensities AVER DENS " << averdens << " -> "
385  << saverdens << " -> " << G4int(1000*averdens) << " " << G4int(1000*averdens)/1000 << " "
386  << G4int(1000*averdens)/1000. << G4endl;
387 
388  G4cout << "GmReadPhantomGeometry::ReadVoxelDensities AVER DENS " << averdens << " -> "
389  << saverdens << " -> " << G4UIcommand::ConvertToString(saverdens) << " Nvoxels "
390  <<(*mppite).second->nvoxels << G4endl;
391  G4String mateName = ((*mppite).first).first->GetName() + "_" +
392  G4UIcommand::ConvertToString(saverdens);
393  fPhantomMaterials.push_back( BuildMaterialWithChangingDensity( (*mppite).first.first,
394  averdens, mateName ) );
395  }
396 
397 }
398 
399 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
400 void DicomPartialDetectorConstruction::ConstructPhantomContainer()
401 {
402  //build a clinder that encompass the partial phantom built in the example
403  // /dicom/intersectWithUserVolume 0. 0. 0. 0.*deg 0. 0. TUBE 0. 200. 100.
404  //---- Extract number of voxels and voxel dimensions
405 
406  G4String contname= "PhantomContainer";
407 
408  //----- Define the volume that contains all the voxels
409  G4Tubs* container_tube = new G4Tubs(contname,0., 200., 100., 0., 360*deg);
410 
412  new G4LogicalVolume( container_tube,
413  fAir,
414  contname,
415  0, 0, 0 );
416 
417  G4ThreeVector posCentreVoxels(0.,0.,0.);
418  //G4cout << " PhantomContainer posCentreVoxels " << posCentreVoxels << G4endl;
420 
422  new G4PVPlacement(rotm, // rotation
423  posCentreVoxels,
424  fContainer_logic, // The logic volume
425  contname, // Name
426  fWorld_logic, // Mother
427  false, // No op. bool.
428  1); // Copy number
429 
430 }
431 
432 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
433 void DicomPartialDetectorConstruction::ConstructPhantom()
434 {
435  G4String OptimAxis = "kXAxis";
436  G4bool bSkipEqualMaterials = 0;
437  G4int regStructureID = 2;
438 
439  G4ThreeVector posCentreVoxels(0.,0.,0.);
440 
441  //----- Build phantom
442  G4String voxelName = "phantom";
443  G4Box* phantom_solid = new G4Box(voxelName, fDimX/2., fDimY/2., fDimZ/2. );
444  G4LogicalVolume* phantom_logic =
445  new G4LogicalVolume( phantom_solid,
446  fAir,
447  voxelName,
448  0, 0, 0 );
449  G4bool pVis = 0;
450  if( !pVis ) {
451  G4VisAttributes* visAtt = new G4VisAttributes( FALSE );
452  phantom_logic->SetVisAttributes( visAtt );
453  }
454 
455  G4double theSmartless = 0.5;
456  // fContainer_logic->SetSmartless( theSmartless );
457  phantom_logic->SetSmartless( theSmartless );
458 
459  fPartialPhantomParam = new G4PartialPhantomParameterisation();
460 
461  fPartialPhantomParam->SetMaterials( fPhantomMaterials );
462  fPartialPhantomParam->SetVoxelDimensions( fDimX/2., fDimY/2., fDimZ/2. );
463  fPartialPhantomParam->SetNoVoxel( fNVoxelX, fNVoxelY, fNVoxelZ );
464  fPartialPhantomParam->SetMaterialIndices( fMateIDs );
465 
466  fPartialPhantomParam->SetFilledIDs(fFilledIDs);
467 
468  fPartialPhantomParam->SetFilledMins(fFilledMins);
469 
470  fPartialPhantomParam->BuildContainerWalls();
471 
472  // G4cout << " Number of Materials " << fPhantomMaterials.size() << G4endl;
473  // G4cout << " SetMaterialIndices(0) " << fMateIDs[0] << G4endl;
474 
475  G4PVParameterised * phantom_phys = 0;
476  if( OptimAxis == "kUndefined" ) {
477  phantom_phys = new G4PVParameterised(voxelName,phantom_logic,fContainer_logic,
478  kUndefined, fNVoxels, fPartialPhantomParam);
479  } else if( OptimAxis == "kXAxis" ) {
480  // G4cout << " optim kX " << G4endl;
481  phantom_phys = new G4PVParameterised(voxelName,phantom_logic,fContainer_logic,
482  // kXAxis, fNVoxelX*fNVoxelY*fNVoxelZ, fPartialPhantomParam);
483  kXAxis, fNVoxels, fPartialPhantomParam);
484  fPartialPhantomParam->SetSkipEqualMaterials(bSkipEqualMaterials);
485  } else {
486  G4Exception("GmReadPhantomGeometry::ConstructPhantom",
487  "",
489  (G4String("Wrong argument to parameter GmReadPhantomGeometry:Phantom:OptimAxis:\
490  Only allowed 'kUndefined' or 'kXAxis', it is: "+OptimAxis).c_str()));
491  }
492 
493  phantom_phys->SetRegularStructureId(regStructureID);
494 
495  // G4cout << " GmReadPhantomGeometry phantom_phys " << phantom_phys << " trans "
496  // << phantom_phys->GetTranslation() << G4endl;
497  //- fPartialPhantomParam->BuildContainerSolid(cont_phys);
498 
499  // G4cout << " GmReadPhantomGeometry::constructPhantom ended " << G4endl;
500 
501 }
fin
Definition: AddMC.C:2
CLHEP::Hep3Vector G4ThreeVector
CLHEP::HepRotation G4RotationMatrix
Definition: G4Box.hh:63
ifstream in
Definition: comparison.C:7
const G4String & GetName() const
Definition: G4Material.hh:176
Definition: G4Tubs.hh:84
G4double GetDensity() const
Definition: G4Material.hh:178
static G4String ConvertToString(G4bool boolVal)
Definition: G4UIcommand.cc:357
int G4int
Definition: G4Types.hh:78
void SetFilledIDs(std::multimap< G4int, G4int > fid)
function g(Y1, Y2, PT2)
Definition: hijing1.383.f:5205
G4Material * BuildMaterialWithChangingDensity(const G4Material *origMate, float density, G4String newMateName)
G4GLOB_DLL std::ostream G4cout
void SetVoxelDimensions(G4double halfx, G4double halfy, G4double halfz)
Definition of the DicomPartialDetectorConstruction class.
bool G4bool
Definition: G4Types.hh:79
G4double iz
Definition: TRTMaterials.hh:39
void SetFilledMins(std::map< G4int, std::map< G4int, G4int > > fmins)
#define FALSE
Definition: globals.hh:52
TString part[npart]
Definition: Style.C:32
static G4double ConvertToDouble(const char *st)
Definition: G4UIcommand.cc:429
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
virtual void SetRegularStructureId(G4int Code)
void SetSmartless(G4double s)
void SetNoVoxel(size_t nx, size_t ny, size_t nz)
G4int first
#define G4endl
Definition: G4ios.hh:61
string fname
Definition: test.py:308
void SetMaterialIndices(size_t *matInd)
double G4double
Definition: G4Types.hh:76
void SetMaterials(std::vector< G4Material * > &mates)
std::vector< G4Material * > fMaterials
std::vector< G4Material * > fOriginalMaterials
#define DBL_MAX
Definition: templates.hh:83
void SetVisAttributes(const G4VisAttributes *pVA)
void SetSkipEqualMaterials(G4bool skip)