ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DicomPartialDetectorConstruction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DicomPartialDetectorConstruction.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
29 //
30 
31 #include "globals.hh"
32 
33 #include "DicomPartialDetectorConstruction.hh"
34 
35 #include "G4Box.hh"
36 #include "G4LogicalVolume.hh"
37 #include "G4VPhysicalVolume.hh"
38 #include "G4PVPlacement.hh"
39 #include "G4PVParameterised.hh"
40 #include "G4Material.hh"
41 #include "G4Element.hh"
42 #include "G4VisAttributes.hh"
43 #include "G4Colour.hh"
44 #include "G4ios.hh"
45 #include "G4PartialPhantomParameterisation.hh"
46 #include "G4NistManager.hh"
47 #include "G4UIcommand.hh"
48 #include "G4Tubs.hh"
49 
50 #include "G4tgbVolumeMgr.hh"
51 #include "G4tgbMaterialMgr.hh"
52 #include "G4tgrMessenger.hh"
53 #include "G4tgrMessenger.hh"
54 
55 #include "G4SystemOfUnits.hh"
56 
57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
58 DicomPartialDetectorConstruction::DicomPartialDetectorConstruction()
59  : DicomDetectorConstruction(),
60  fPartialPhantomParam(0),
61  fNVoxels(0),
62  fDimX(0),
63  fDimY(0),
64  fDimZ(0),
65  fOffsetX(0),
66  fOffsetY(0),
67  fOffsetZ(0)
68 {
69 }
70 
71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
72 DicomPartialDetectorConstruction::~DicomPartialDetectorConstruction()
73 {
74 }
75 
76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
77 G4VPhysicalVolume* DicomPartialDetectorConstruction::Construct()
78 {
79  InitialisationOfMaterials();
80 
81  //----- Build world
82  G4double worldXDimension = 1.*m;
83  G4double worldYDimension = 1.*m;
84  G4double worldZDimension = 1.*m;
85 
86  G4Box* world_box = new G4Box( "WorldSolid",
87  worldXDimension,
88  worldYDimension,
89  worldZDimension );
90 
91  fWorld_logic = new G4LogicalVolume( world_box,
92  fAir,
93  "WorldLogical",
94  0, 0, 0 );
95 
96  G4VPhysicalVolume* world_phys = new G4PVPlacement( 0,
97  G4ThreeVector(0,0,0),
98  "World",
99  fWorld_logic,
100  0,
101  false,
102  0 );
103 
104  ReadPhantomData();
105 
106  ConstructPhantomContainer();
107  ConstructPhantom();
108 
109  return world_phys;
110 }
111 
112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
113 void DicomPartialDetectorConstruction::ReadPhantomData()
114 {
115  G4String dataFile = "Data.dat";
116  std::ifstream finDF(dataFile.c_str());
117  G4String fname;
118  if(finDF.good() != 1 ) {
119  G4Exception(" DicomPartialDetectorConstruction::ReadPhantomData",
120  "",
121  FatalErrorInArgument,
122  " Problem reading data file: Data.dat");
123  }
124 
125  G4int compression;
126  finDF >> compression; // not used here
127 
128  G4int numFiles;
129  finDF >> numFiles; // only 1 file supported for the moment
130  if( numFiles != 1 ) {
131  G4Exception("DicomPartialDetectorConstruction::ReadPhantomData",
132  "",
133  FatalErrorInArgument,
134  "More than 1 DICOM file is not supported");
135  }
136  for(G4int i = 0; i < numFiles; i++ ) {
137  finDF >> fname;
138  //--- Read one data file
139  fname += ".g4pdcm";
140  ReadPhantomDataFile(fname);
141  }
142 
143  finDF.close();
144 }
145 
146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
147 void DicomPartialDetectorConstruction::
148 ReadPhantomDataFile(const G4String& fname)
149 {
150  // G4String filename = "phantom.g4pdcm";
151 #ifdef G4VERBOSE
152  G4cout << " DicomDetectorConstruction::ReadPhantomDataFile opening file "
153  << fname << G4endl;
154 #endif
155  std::ifstream fin(fname.c_str(), std::ios_base::in);
156  if( !fin.is_open() ) {
157  G4Exception("DicomPartialDetectorConstruction::ReadPhantomDataFile",
158  "",
159  FatalException,
160  G4String("File not found " + fname).c_str());
161  }
162  G4int nMaterials;
163  fin >> nMaterials;
164  G4String stemp;
165  G4int nmate;
166  for( G4int ii = 0; ii < nMaterials; ii++ ){
167  fin >> nmate >> stemp;
168  G4cout << "DicomPartialDetectorConstruction::ReadPhantomData reading nmate "
169  << ii << " = " << nmate << " mate " << stemp << G4endl;
170  if( ii != nmate )
171  G4Exception("DicomPartialDetectorConstruction::ReadPhantomData",
172  "",
173  FatalErrorInArgument,
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  G4int 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::pair<G4int,G4int>(ifxdiff+fNVoxels-1, ifxmin1));
212  G4cout << "DicomPartialDetectorConstruction::ReadPhantomData insert "
213  << " FilledIDs " << ifxdiff+fNVoxels-1 << " min " << ifxmin1
214  << " N= " << fFilledIDs.size() << G4endl;
215  //filledIDs[iz*fNVoxelY+iy+1] = ifxmax1-ifxmin1+1;
216  for( G4int ix = 0; ix < fNVoxelX; ix++ ) {
217  if( ix >= G4int(ifxmin1) && ix <= G4int(ifxmax1) ) {
218  fNVoxels++;
219  } else {
220  }
221  }
222  }
223  fFilledMins[iz] = ifmin;
224  fFilledMaxs[iz] = ifmax;
225  }
226 
227  //--- Read material IDs
228  G4int mateID1;
229  fMateIDs = new size_t[fNVoxelX*fNVoxelY*fNVoxelZ];
230  G4int copyNo = 0;
231  for( G4int iz = 0; iz < fNVoxelZ; iz++ ) {
232  std::map< G4int, G4int > ifmin = fFilledMins[iz];
233  std::map< G4int, G4int > ifmax = fFilledMaxs[iz];
234  for( G4int iy = 0; iy < fNVoxelY; iy++ ) {
235  ifxmin1 = ifmin[iy];
236  ifxmax1 = ifmax[iy];
237  for( G4int ix = 0; ix < fNVoxelX; ix++ ) {
238  if( ix >= G4int(ifxmin1) && ix <= G4int(ifxmax1) ) {
239  fin >> mateID1;
240  if( mateID1 < 0 || mateID1 >= nMaterials ) {
241  G4Exception("DicomPartialDetectorConstruction::ReadPhantomData",
242  "",
243  FatalException,
244  G4String("Wrong index in phantom file: It should be between 0 and "
245  + G4UIcommand::ConvertToString(nMaterials-1)
246  + ", while it is "
247  + G4UIcommand::ConvertToString(mateID1)).c_str());
248  }
249  fMateIDs[copyNo] = mateID1;
250  copyNo++;
251  }
252  }
253  }
254  }
255 
256  ReadVoxelDensitiesPartial( fin );
257 
258  fin.close();
259 }
260 
261 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
262 void DicomPartialDetectorConstruction::
263 ReadVoxelDensitiesPartial( std::ifstream& fin )
264 {
265  std::map<G4int, std::pair<G4double,G4double> > densiMinMax;
266  std::map<G4int, std::pair<G4double,G4double> >::iterator mpite;
267  for( G4int ii = 0; ii < G4int(fOriginalMaterials.size()); ++ii )
268  {
269  densiMinMax[ii] = std::pair<G4double,G4double>(DBL_MAX,-DBL_MAX);
270  }
271 
272  //----- Define density differences (maximum density difference to create
273  // a new material)
274  char* part = std::getenv( "DICOM_CHANGE_MATERIAL_DENSITY" );
275  G4double densityDiff = -1.;
276  if( part ) densityDiff = G4UIcommand::ConvertToDouble(part);
277  std::map<G4int,G4double> densityDiffs;
278  if( densityDiff != -1. )
279  {
280  for( G4int ii = 0; ii < G4int(fOriginalMaterials.size()); ++ii )
281  {
282  densityDiffs[ii] = densityDiff; //currently all materials
283  //with same step
284  }
285  }
286  else
287  {
288  if( fMaterials.size() == 0 ) { // do it only for first slice
289  for( unsigned int ii = 0; ii < fOriginalMaterials.size(); ++ii )
290  {
291  fMaterials.push_back( fOriginalMaterials[ii] );
292  }
293  }
294  }
295  // densityDiffs[0] = 0.0001; //fAir
296 
297  //--- Calculate the average material density for each material/density bin
298  std::map< std::pair<G4Material*,G4int>, matInfo* > newMateDens;
299 
300  G4double dens1;
301  G4int ifxmin1, ifxmax1;
302 
303  //---- Read the material densities
304  G4int copyNo = 0;
305  for( G4int iz = 0; iz < fNVoxelZ; ++iz )
306  {
307  std::map< G4int, G4int > ifmin = fFilledMins[iz];
308  std::map< G4int, G4int > ifmax = fFilledMaxs[iz];
309  for( G4int iy = 0; iy < fNVoxelY; ++iy )
310  {
311  ifxmin1 = ifmin[iy];
312  ifxmax1 = ifmax[iy];
313  for( G4int ix = 0; ix < fNVoxelX; ++ix )
314  {
315  if( ix >= G4int(ifxmin1) && ix <= G4int(ifxmax1) ) {
316  fin >> dens1;
317  //--- store the minimum and maximum density for each material
318  //(just for printing)
319  mpite = densiMinMax.find( G4int(fMateIDs[copyNo]) );
320  if( dens1 < (*mpite).second.first ) (*mpite).second.first = dens1;
321  if( dens1 > (*mpite).second.second ) (*mpite).second.second = dens1;
322 
323  //--- Get material from original list of material in file
324  G4int mateID = G4int(fMateIDs[copyNo]);
325  G4Material* mate = fOriginalMaterials[mateID];
326  // G4cout << copyNo << " mateID " << mateID << G4endl;
327  //--- Check if density is equal to the original material density
328  if( std::fabs(dens1 - mate->GetDensity()/g*cm3 ) < 1.e-9 ) continue;
329 
330  //--- Build material name with fOriginalMaterials name + density
331  // float densityBin = densityDiffs[mateID] *
332  // (G4int(dens1/densityDiffs[mateID])+0.5);
333  G4String newMateName = mate->GetName();
334 
335  G4int densityBin = 0;
336  // G4cout << " densityBin " << densityBin << " "
337  // << dens1 << " " <<densityDiffs[mateID] << G4endl;
338 
339  if( densityDiff != -1.) {
340  densityBin = (G4int(dens1/densityDiffs[mateID]));
341  newMateName =
342  mate->GetName()+G4UIcommand::ConvertToString(densityBin);
343  }
344 
345  //--- Look if it is the first voxel with this material/densityBin
346  std::pair<G4Material*,G4int> matdens(mate, densityBin );
347 
348  auto mppite = newMateDens.find( matdens );
349  if( mppite != newMateDens.cend() ){
350  matInfo* mi = (*mppite).second;
351  mi->fSumdens += dens1;
352  mi->fNvoxels++;
353  fMateIDs[copyNo] = fOriginalMaterials.size()-1 + mi->fId;
354  // G4cout << copyNo << " mat new again "
355  //<< fOriginalMaterials.size()-1 + mi->id << " " << mi->id << G4endl;
356  } else {
357  matInfo* mi = new matInfo;
358  mi->fSumdens = dens1;
359  mi->fNvoxels = 1;
360  mi->fId = G4int(newMateDens.size()+1);
361  newMateDens[matdens] = mi;
362  fMateIDs[copyNo] = fOriginalMaterials.size()-1 + mi->fId;
363  // G4cout << copyNo << " mat new first "
364  // << fOriginalMaterials.size()-1 + mi->fId << G4endl;
365  }
366  copyNo++;
367  // G4cout << ix << " " << iy << " " << iz
368  // << " filling fMateIDs " << copyNo << " = " << atoi(cid)-1 << G4endl;
369  // fMateIDs[copyNo] = atoi(cid)-1;
370  }
371  }
372  }
373  }
374 
375  //----- Build the list of phantom materials that go to Parameterisation
376  //--- Add original materials
377  for( auto mimite = fOriginalMaterials.cbegin();
378  mimite != fOriginalMaterials.cend(); ++mimite )
379  {
380  fPhantomMaterials.push_back( (*mimite) );
381  }
382  //
383  //---- Build and add new materials
384  for( auto mppite= newMateDens.cbegin();
385  mppite != newMateDens.cend(); ++mppite )
386  {
387  G4double averdens = (*mppite).second->fSumdens/(*mppite).second->fNvoxels;
388  G4double saverdens = G4int(1000.001*averdens)/1000.;
389  G4cout << "GmReadPhantomGeometry::ReadVoxelDensities AVER DENS "
390  << averdens << " -> " << saverdens << " -> "
391  << G4int(1000*averdens) << " " << G4int(1000*averdens)/1000 << " "
392  << G4int(1000*averdens)/1000. << G4endl;
393 
394  G4cout << "GmReadPhantomGeometry::ReadVoxelDensities AVER DENS "
395  << averdens << " -> " << saverdens << " -> "
396  << G4UIcommand::ConvertToString(saverdens) << " Nvoxels "
397  <<(*mppite).second->fNvoxels << G4endl;
398  G4String mateName = ((*mppite).first).first->GetName() + "_" +
399  G4UIcommand::ConvertToString(saverdens);
400  fPhantomMaterials.push_back( BuildMaterialWithChangingDensity(
401  (*mppite).first.first,
402  G4float(averdens), mateName ) );
403  }
404 }
405 
406 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
407 void DicomPartialDetectorConstruction::ConstructPhantomContainer()
408 {
409  //build a clinder that encompass the partial phantom built in the example
410  // /dicom/intersectWithUserVolume 0. 0. 0. 0.*deg 0. 0. TUBE 0. 200. 100.
411  //---- Extract number of voxels and voxel dimensions
412 
413  G4String contname= "PhantomContainer";
414 
415  //----- Define the volume that contains all the voxels
416  G4Tubs* container_tube = new G4Tubs(contname,0., 200., 100., 0., 360*deg);
417 
418  fContainer_logic =
419  new G4LogicalVolume( container_tube,
420  fAir,
421  contname,
422  0, 0, 0 );
423 
424  G4ThreeVector posCentreVoxels(0.,0.,0.);
425  //G4cout << " PhantomContainer posCentreVoxels " << posCentreVoxels << G4endl;
426  G4RotationMatrix* rotm = new G4RotationMatrix;
427 
428  fContainer_phys =
429  new G4PVPlacement(rotm, // rotation
430  posCentreVoxels,
431  fContainer_logic, // The logic volume
432  contname, // Name
433  fWorld_logic, // Mother
434  false, // No op. bool.
435  1); // Copy number
436 
437 }
438 
439 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
440 void DicomPartialDetectorConstruction::ConstructPhantom()
441 {
442  G4String OptimAxis = "kXAxis";
443  G4bool bSkipEqualMaterials = 0;
444  G4int regStructureID = 2;
445 
446  G4ThreeVector posCentreVoxels(0.,0.,0.);
447 
448  //----- Build phantom
449  G4String voxelName = "phantom";
450  G4Box* phantom_solid = new G4Box(voxelName, fDimX/2., fDimY/2., fDimZ/2. );
451  G4LogicalVolume* phantom_logic =
452  new G4LogicalVolume( phantom_solid,
453  fAir,
454  voxelName,
455  0, 0, 0 );
456  G4bool pVis = 0;
457  if( !pVis ) {
458  G4VisAttributes* visAtt = new G4VisAttributes( FALSE );
459  phantom_logic->SetVisAttributes( visAtt );
460  }
461 
462  G4double theSmartless = 0.5;
463  // fContainer_logic->SetSmartless( theSmartless );
464  phantom_logic->SetSmartless( theSmartless );
465 
466  fPartialPhantomParam = new G4PartialPhantomParameterisation();
467 
468  fPartialPhantomParam->SetMaterials( fPhantomMaterials );
469  fPartialPhantomParam->SetVoxelDimensions( fDimX/2., fDimY/2., fDimZ/2. );
470  fPartialPhantomParam->SetNoVoxel( fNVoxelX, fNVoxelY, fNVoxelZ );
471  fPartialPhantomParam->SetMaterialIndices( fMateIDs );
472 
473  fPartialPhantomParam->SetFilledIDs(fFilledIDs);
474 
475  fPartialPhantomParam->SetFilledMins(fFilledMins);
476 
477  fPartialPhantomParam->BuildContainerWalls();
478 
479  // G4cout << " Number of Materials " << fPhantomMaterials.size() << G4endl;
480  // G4cout << " SetMaterialIndices(0) " << fMateIDs[0] << G4endl;
481 
482  G4PVParameterised * phantom_phys = 0;
483  if( OptimAxis == "kUndefined" ) {
484  phantom_phys =
485  new G4PVParameterised(voxelName,phantom_logic,fContainer_logic,
486  kUndefined, fNVoxels, fPartialPhantomParam);
487  } else if( OptimAxis == "kXAxis" ) {
488  // G4cout << " optim kX " << G4endl;
489  phantom_phys =
490  new G4PVParameterised(voxelName,phantom_logic,fContainer_logic,
491  kXAxis, fNVoxels, fPartialPhantomParam);
492  fPartialPhantomParam->SetSkipEqualMaterials(bSkipEqualMaterials);
493  } else {
494  G4Exception("GmReadPhantomGeometry::ConstructPhantom",
495  "",
496  FatalErrorInArgument,
497  G4String("Wrong argument to parameter \
498  GmReadPhantomGeometry:Phantom:OptimAxis: \
499  Only allowed 'kUndefined' or 'kXAxis', it is: "+OptimAxis).c_str());
500  }
501 
502  phantom_phys->SetRegularStructureId(regStructureID);
503 
504 }