Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/medical/DICOM/src/DicomDetectorConstruction.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /examples/extended/medical/DICOM/src/DicomDetectorConstruction.cc (Version 11.3.0) and /examples/extended/medical/DICOM/src/DicomDetectorConstruction.cc (Version 8.0)


  1 //                                                  1 
  2 // *******************************************    
  3 // * License and Disclaimer                       
  4 // *                                              
  5 // * The  Geant4 software  is  copyright of th    
  6 // * the Geant4 Collaboration.  It is provided    
  7 // * conditions of the Geant4 Software License    
  8 // * LICENSE and available at  http://cern.ch/    
  9 // * include a list of copyright holders.         
 10 // *                                              
 11 // * Neither the authors of this software syst    
 12 // * institutes,nor the agencies providing fin    
 13 // * work  make  any representation or  warran    
 14 // * regarding  this  software system or assum    
 15 // * use.  Please see the license in the file     
 16 // * for the full disclaimer and the limitatio    
 17 // *                                              
 18 // * This  code  implementation is the result     
 19 // * technical work of the GEANT4 collaboratio    
 20 // * By using,  copying,  modifying or  distri    
 21 // * any work based  on the software)  you  ag    
 22 // * use  in  resulting  scientific  publicati    
 23 // * acceptance of all terms of the Geant4 Sof    
 24 // *******************************************    
 25 //                                                
 26 //                                                
 27 /// \file  medical/DICOM/src/DicomDetectorCons    
 28 /// \brief Implementation of the DicomDetector    
 29 //                                                
 30                                                   
 31 #include "DicomDetectorConstruction.hh"           
 32                                                   
 33 #include "CLHEP/Units/SystemOfUnits.h"            
 34 #include "DicomHandler.hh"                        
 35 #include "DicomPhantomZSliceHeader.hh"            
 36                                                   
 37 #include "G4Box.hh"                               
 38 #include "G4Element.hh"                           
 39 #include "G4LogicalVolume.hh"                     
 40 #include "G4Material.hh"                          
 41 #include "G4NistManager.hh"                       
 42 #include "G4PVPlacement.hh"                       
 43 #include "G4PhysicalConstants.hh"                 
 44 #include "G4UIcommand.hh"                         
 45 #include "G4VPhysicalVolume.hh"                   
 46 #include "globals.hh"                             
 47                                                   
 48 #ifdef G4_DCMTK                                   
 49 #  include "DicomFileMgr.hh"                      
 50 #endif                                            
 51 #include "G4VisAttributes.hh"                     
 52                                                   
 53 using CLHEP::cm3;                                 
 54 using CLHEP::g;                                   
 55 using CLHEP::m;                                   
 56 using CLHEP::mg;                                  
 57 using CLHEP::mole;                                
 58 using CLHEP::perCent;                             
 59                                                   
 60 //....oooOO0OOooo........oooOO0OOooo........oo    
 61 DicomDetectorConstruction::DicomDetectorConstr    
 62   : G4VUserDetectorConstruction(),                
 63     fAir(0),                                      
 64                                                   
 65     fWorld_solid(0),                              
 66     fWorld_logic(0),                              
 67     fWorld_phys(0),                               
 68                                                   
 69     fContainer_solid(0),                          
 70     fContainer_logic(0),                          
 71     fContainer_phys(0),                           
 72                                                   
 73     fNoFiles(0),                                  
 74     fMateIDs(0),                                  
 75                                                   
 76     fZSliceHeaderMerged(0),                       
 77                                                   
 78     fNoVoxelsX(0),                                
 79     fNoVoxelsY(0),                                
 80     fNoVoxelsZ(0),                                
 81     fVoxelHalfDimX(0),                            
 82     fVoxelHalfDimY(0),                            
 83     fVoxelHalfDimZ(0),                            
 84                                                   
 85     fConstructed(false)                           
 86 {}                                                
 87                                                   
 88 //....oooOO0OOooo........oooOO0OOooo........oo    
 89 DicomDetectorConstruction::~DicomDetectorConst    
 90                                                   
 91 //....oooOO0OOooo........oooOO0OOooo........oo    
 92 G4VPhysicalVolume* DicomDetectorConstruction::    
 93 {                                                 
 94   if (!fConstructed || fWorld_phys == 0) {        
 95     fConstructed = true;                          
 96     InitialisationOfMaterials();                  
 97                                                   
 98     //----- Build world                           
 99     G4double worldXDimension = 1. * m;            
100     G4double worldYDimension = 1. * m;            
101     G4double worldZDimension = 1. * m;            
102                                                   
103     fWorld_solid = new G4Box("WorldSolid", wor    
104                                                   
105     fWorld_logic = new G4LogicalVolume(fWorld_    
106                                                   
107     fWorld_phys = new G4PVPlacement(0, G4Three    
108                                                   
109     fWorld_logic->SetVisAttributes(G4VisAttrib    
110                                                   
111 #ifdef G4_DCMTK                                   
112     ReadPhantomDataNew();                         
113     ConstructPhantomContainerNew();               
114 #else                                             
115     ReadPhantomData();                            
116     ConstructPhantomContainer();                  
117 #endif                                            
118                                                   
119     ConstructPhantom();                           
120   }                                               
121   return fWorld_phys;                             
122 }                                                 
123                                                   
124 //....oooOO0OOooo........oooOO0OOooo........oo    
125 void DicomDetectorConstruction::Initialisation    
126 {                                                 
127   // Creating elements :                          
128   G4double z, a, density;                         
129   G4String name, symbol;                          
130                                                   
131   G4Element* elC = new G4Element(name = "Carbo    
132   G4Element* elH = new G4Element(name = "Hydro    
133   G4Element* elN = new G4Element(name = "Nitro    
134   G4Element* elO = new G4Element(name = "Oxyge    
135   G4Element* elNa =                               
136     new G4Element(name = "Sodium", symbol = "N    
137   G4Element* elMg =                               
138     new G4Element(name = "Magnesium", symbol =    
139   G4Element* elP =                                
140     new G4Element(name = "Phosphorus", symbol     
141   G4Element* elS = new G4Element(name = "Sulfu    
142   G4Element* elCl =                               
143     new G4Element(name = "Chlorine", symbol =     
144   G4Element* elK =                                
145     new G4Element(name = "Potassium", symbol =    
146                                                   
147   G4Element* elFe = new G4Element(name = "Iron    
148                                                   
149   G4Element* elCa = new G4Element(name = "Calc    
150                                                   
151   G4Element* elZn = new G4Element(name = "Zinc    
152                                                   
153   // Creating Materials :                         
154   G4int numberofElements;                         
155                                                   
156   // Air                                          
157   fAir = new G4Material("Air", 1.290 * mg / cm    
158   fAir->AddElement(elN, 0.7);                     
159   fAir->AddElement(elO, 0.3);                     
160                                                   
161   // Soft tissue (ICRP - NIST)                    
162   G4Material* softTissue = new G4Material("Sof    
163   softTissue->AddElement(elH, 10.4472 * perCen    
164   softTissue->AddElement(elC, 23.219 * perCent    
165   softTissue->AddElement(elN, 2.488 * perCent)    
166   softTissue->AddElement(elO, 63.0238 * perCen    
167   softTissue->AddElement(elNa, 0.113 * perCent    
168   softTissue->AddElement(elMg, 0.0113 * perCen    
169   softTissue->AddElement(elP, 0.113 * perCent)    
170   softTissue->AddElement(elS, 0.199 * perCent)    
171   softTissue->AddElement(elCl, 0.134 * perCent    
172   softTissue->AddElement(elK, 0.199 * perCent)    
173   softTissue->AddElement(elCa, 0.023 * perCent    
174   softTissue->AddElement(elFe, 0.005 * perCent    
175   softTissue->AddElement(elZn, 0.003 * perCent    
176                                                   
177   //  Lung Inhale                                 
178   G4Material* lunginhale =                        
179     new G4Material("LungInhale", density = 0.2    
180   lunginhale->AddElement(elH, 0.103);             
181   lunginhale->AddElement(elC, 0.105);             
182   lunginhale->AddElement(elN, 0.031);             
183   lunginhale->AddElement(elO, 0.749);             
184   lunginhale->AddElement(elNa, 0.002);            
185   lunginhale->AddElement(elP, 0.002);             
186   lunginhale->AddElement(elS, 0.003);             
187   lunginhale->AddElement(elCl, 0.002);            
188   lunginhale->AddElement(elK, 0.003);             
189                                                   
190   // Lung exhale                                  
191   G4Material* lungexhale =                        
192     new G4Material("LungExhale", density = 0.5    
193   lungexhale->AddElement(elH, 0.103);             
194   lungexhale->AddElement(elC, 0.105);             
195   lungexhale->AddElement(elN, 0.031);             
196   lungexhale->AddElement(elO, 0.749);             
197   lungexhale->AddElement(elNa, 0.002);            
198   lungexhale->AddElement(elP, 0.002);             
199   lungexhale->AddElement(elS, 0.003);             
200   lungexhale->AddElement(elCl, 0.002);            
201   lungexhale->AddElement(elK, 0.003);             
202                                                   
203   // Adipose tissue                               
204   G4Material* adiposeTissue =                     
205     new G4Material("AdiposeTissue", density =     
206   adiposeTissue->AddElement(elH, 0.114);          
207   adiposeTissue->AddElement(elC, 0.598);          
208   adiposeTissue->AddElement(elN, 0.007);          
209   adiposeTissue->AddElement(elO, 0.278);          
210   adiposeTissue->AddElement(elNa, 0.001);         
211   adiposeTissue->AddElement(elS, 0.001);          
212   adiposeTissue->AddElement(elCl, 0.001);         
213                                                   
214   // Brain (ICRP - NIST)                          
215   G4Material* brainTissue = new G4Material("Br    
216   brainTissue->AddElement(elH, 11.0667 * perCe    
217   brainTissue->AddElement(elC, 12.542 * perCen    
218   brainTissue->AddElement(elN, 1.328 * perCent    
219   brainTissue->AddElement(elO, 73.7723 * perCe    
220   brainTissue->AddElement(elNa, 0.1840 * perCe    
221   brainTissue->AddElement(elMg, 0.015 * perCen    
222   brainTissue->AddElement(elP, 0.356 * perCent    
223   brainTissue->AddElement(elS, 0.177 * perCent    
224   brainTissue->AddElement(elCl, 0.236 * perCen    
225   brainTissue->AddElement(elK, 0.31 * perCent)    
226   brainTissue->AddElement(elCa, 0.009 * perCen    
227   brainTissue->AddElement(elFe, 0.005 * perCen    
228   brainTissue->AddElement(elZn, 0.001 * perCen    
229                                                   
230   // Breast                                       
231   G4Material* breast = new G4Material("Breast"    
232   breast->AddElement(elH, 0.109);                 
233   breast->AddElement(elC, 0.506);                 
234   breast->AddElement(elN, 0.023);                 
235   breast->AddElement(elO, 0.358);                 
236   breast->AddElement(elNa, 0.001);                
237   breast->AddElement(elP, 0.001);                 
238   breast->AddElement(elS, 0.001);                 
239   breast->AddElement(elCl, 0.001);                
240                                                   
241   // Spinal Disc                                  
242   G4Material* spinalDisc = new G4Material("Spi    
243   spinalDisc->AddElement(elH, 9.60 * perCent);    
244   spinalDisc->AddElement(elC, 9.90 * perCent);    
245   spinalDisc->AddElement(elN, 2.20 * perCent);    
246   spinalDisc->AddElement(elO, 74.40 * perCent)    
247   spinalDisc->AddElement(elNa, 0.50 * perCent)    
248   spinalDisc->AddElement(elP, 2.20 * perCent);    
249   spinalDisc->AddElement(elS, 0.90 * perCent);    
250   spinalDisc->AddElement(elCl, 0.30 * perCent)    
251                                                   
252   // Water                                        
253   G4Material* water = new G4Material("Water",     
254   water->AddElement(elH, 0.112);                  
255   water->AddElement(elO, 0.888);                  
256                                                   
257   // Muscle                                       
258   G4Material* muscle = new G4Material("Muscle"    
259   muscle->AddElement(elH, 0.102);                 
260   muscle->AddElement(elC, 0.143);                 
261   muscle->AddElement(elN, 0.034);                 
262   muscle->AddElement(elO, 0.710);                 
263   muscle->AddElement(elNa, 0.001);                
264   muscle->AddElement(elP, 0.002);                 
265   muscle->AddElement(elS, 0.003);                 
266   muscle->AddElement(elCl, 0.001);                
267   muscle->AddElement(elK, 0.004);                 
268                                                   
269   // Liver                                        
270   G4Material* liver = new G4Material("Liver",     
271   liver->AddElement(elH, 0.102);                  
272   liver->AddElement(elC, 0.139);                  
273   liver->AddElement(elN, 0.030);                  
274   liver->AddElement(elO, 0.716);                  
275   liver->AddElement(elNa, 0.002);                 
276   liver->AddElement(elP, 0.003);                  
277   liver->AddElement(elS, 0.003);                  
278   liver->AddElement(elCl, 0.002);                 
279   liver->AddElement(elK, 0.003);                  
280                                                   
281   // Tooth Dentin                                 
282   G4Material* toothDentin = new G4Material("To    
283   toothDentin->AddElement(elH, 2.67 * perCent)    
284   toothDentin->AddElement(elC, 12.77 * perCent    
285   toothDentin->AddElement(elN, 4.27 * perCent)    
286   toothDentin->AddElement(elO, 40.40 * perCent    
287   toothDentin->AddElement(elNa, 0.65 * perCent    
288   toothDentin->AddElement(elMg, 0.59 * perCent    
289   toothDentin->AddElement(elP, 11.86 * perCent    
290   toothDentin->AddElement(elCl, 0.04 * perCent    
291   toothDentin->AddElement(elCa, 26.74 * perCen    
292   toothDentin->AddElement(elZn, 0.01 * perCent    
293                                                   
294   // Trabecular Bone                              
295   G4Material* trabecularBone =                    
296     new G4Material("TrabecularBone", density =    
297   trabecularBone->AddElement(elH, 0.085);         
298   trabecularBone->AddElement(elC, 0.404);         
299   trabecularBone->AddElement(elN, 0.058);         
300   trabecularBone->AddElement(elO, 0.367);         
301   trabecularBone->AddElement(elNa, 0.001);        
302   trabecularBone->AddElement(elMg, 0.001);        
303   trabecularBone->AddElement(elP, 0.034);         
304   trabecularBone->AddElement(elS, 0.002);         
305   trabecularBone->AddElement(elCl, 0.002);        
306   trabecularBone->AddElement(elK, 0.001);         
307   trabecularBone->AddElement(elCa, 0.044);        
308   trabecularBone->AddElement(elFe, 0.001);        
309                                                   
310   // Trabecular bone used in the DICOM Head       
311                                                   
312   G4Material* trabecularBone_head =               
313     new G4Material("TrabecularBone_HEAD", 1.18    
314   trabecularBone_head->AddElement(elH, 8.50 *     
315   trabecularBone_head->AddElement(elC, 40.40 *    
316   trabecularBone_head->AddElement(elN, 2.80 *     
317   trabecularBone_head->AddElement(elO, 36.70 *    
318   trabecularBone_head->AddElement(elNa, 0.10 *    
319   trabecularBone_head->AddElement(elMg, 0.10 *    
320   trabecularBone_head->AddElement(elP, 3.40 *     
321   trabecularBone_head->AddElement(elS, 0.20 *     
322   trabecularBone_head->AddElement(elCl, 0.20 *    
323   trabecularBone_head->AddElement(elK, 0.10 *     
324   trabecularBone_head->AddElement(elCa, 7.40 *    
325   trabecularBone_head->AddElement(elFe, 0.10 *    
326                                                   
327   // Dense Bone                                   
328   G4Material* denseBone =                         
329     new G4Material("DenseBone", density = 1.57    
330   denseBone->AddElement(elH, 0.056);              
331   denseBone->AddElement(elC, 0.235);              
332   denseBone->AddElement(elN, 0.050);              
333   denseBone->AddElement(elO, 0.434);              
334   denseBone->AddElement(elNa, 0.001);             
335   denseBone->AddElement(elMg, 0.001);             
336   denseBone->AddElement(elP, 0.072);              
337   denseBone->AddElement(elS, 0.003);              
338   denseBone->AddElement(elCl, 0.001);             
339   denseBone->AddElement(elK, 0.001);              
340   denseBone->AddElement(elCa, 0.146);             
341                                                   
342   // Cortical Bone (ICRP - NIST)                  
343   G4Material* corticalBone = new G4Material("C    
344   corticalBone->AddElement(elH, 4.7234 * perCe    
345   corticalBone->AddElement(elC, 14.4330 * perC    
346   corticalBone->AddElement(elN, 4.199 * perCen    
347   corticalBone->AddElement(elO, 44.6096 * perC    
348   corticalBone->AddElement(elMg, 0.22 * perCen    
349   corticalBone->AddElement(elP, 10.497 * perCe    
350   corticalBone->AddElement(elS, 0.315 * perCen    
351   corticalBone->AddElement(elCa, 20.993 * perC    
352   corticalBone->AddElement(elZn, 0.01 * perCen    
353                                                   
354   // Tooth enamel                                 
355   G4Material* toothEnamel = new G4Material("To    
356   toothEnamel->AddElement(elH, 0.95 * perCent)    
357   toothEnamel->AddElement(elC, 1.11 * perCent)    
358   toothEnamel->AddElement(elN, 0.23 * perCent)    
359   toothEnamel->AddElement(elO, 41.66 * perCent    
360   toothEnamel->AddElement(elNa, 0.79 * perCent    
361   toothEnamel->AddElement(elMg, 0.23 * perCent    
362   toothEnamel->AddElement(elP, 18.71 * perCent    
363   toothEnamel->AddElement(elCl, 0.34 * perCent    
364   toothEnamel->AddElement(elCa, 35.97 * perCen    
365   toothEnamel->AddElement(elZn, 0.02 * perCent    
366                                                   
367 #ifdef DICOM_USE_HEAD                             
368   //----- Put the materials in a vector HEAD P    
369   fOriginalMaterials.push_back(fAir);  // 0.00    
370   fOriginalMaterials.push_back(softTissue);  /    
371   fOriginalMaterials.push_back(brainTissue);      
372   fOriginalMaterials.push_back(spinalDisc);  /    
373   fOriginalMaterials.push_back(trabecularBone_    
374   fOriginalMaterials.push_back(toothDentin);      
375   fOriginalMaterials.push_back(corticalBone);     
376   fOriginalMaterials.push_back(toothEnamel);      
377   G4cout << "The materials of the DICOM Head h    
378 #else                                             
379   fOriginalMaterials.push_back(fAir);  // rho     
380   fOriginalMaterials.push_back(lunginhale);  /    
381   fOriginalMaterials.push_back(lungexhale);  /    
382   fOriginalMaterials.push_back(adiposeTissue);    
383   fOriginalMaterials.push_back(breast);  // rh    
384   fOriginalMaterials.push_back(water);  // rho    
385   fOriginalMaterials.push_back(muscle);  // rh    
386   fOriginalMaterials.push_back(liver);  // rho    
387   fOriginalMaterials.push_back(trabecularBone)    
388   fOriginalMaterials.push_back(denseBone);  //    
389   G4cout << "Default materials of the DICOM Ex    
390 #endif                                            
391 }                                                 
392                                                   
393 //....oooOO0OOooo........oooOO0OOooo........oo    
394 void DicomDetectorConstruction::ReadPhantomDat    
395 {                                                 
396 #ifdef G4_DCMTK                                   
397   G4String fileName = DicomFileMgr::GetInstanc    
398                                                   
399   std::ifstream fin(fileName);                    
400   std::vector<G4String> wl;                       
401   G4int nMaterials;                               
402   fin >> nMaterials;                              
403   G4String mateName;                              
404   G4int nmate;                                    
405   for (G4int ii = 0; ii < nMaterials; ++ii) {     
406     fin >> nmate;                                 
407     fin >> mateName;                              
408     if (mateName[0] == '"' && mateName[mateNam    
409       mateName = mateName.substr(1, mateName.l    
410     }                                             
411     G4cout << "GmReadPhantomG4Geometry::ReadPh    
412            << " mate " << mateName << G4endl;     
413     if (ii != nmate)                              
414       G4Exception("GmReadPhantomG4Geometry::Re    
415                   FatalErrorInArgument,           
416                   "Material number should be i    
417                                                   
418     G4Material* mate = 0;                         
419     const G4MaterialTable* matTab = G4Material    
420     for (auto matite = matTab->cbegin(); matit    
421       if ((*matite)->GetName() == mateName) {     
422         mate = *matite;                           
423       }                                           
424     }                                             
425     if (mate == 0) {                              
426       mate = G4NistManager::Instance()->FindOr    
427     }                                             
428     if (!mate)                                    
429       G4Exception("GmReadPhantomG4Geometry::Re    
430                   FatalErrorInArgument, ("Mate    
431     fPhantomMaterialsOriginal[nmate] = mate;      
432   }                                               
433                                                   
434   fin >> fNoVoxelsX >> fNoVoxelsY >> fNoVoxels    
435   G4cout << "GmReadPhantomG4Geometry::ReadPhan    
436          << fNoVoxelsY << " " << fNoVoxelsZ <<    
437   fin >> fMinX >> fMaxX;                          
438   fin >> fMinY >> fMaxY;                          
439   fin >> fMinZ >> fMaxZ;                          
440   fVoxelHalfDimX = (fMaxX - fMinX) / fNoVoxels    
441   fVoxelHalfDimY = (fMaxY - fMinY) / fNoVoxels    
442   fVoxelHalfDimZ = (fMaxZ - fMinZ) / fNoVoxels    
443 #  ifdef G4VERBOSE                                
444   G4cout << " Extension in X " << fMinX << " "    
445          << " " << fMaxY << G4endl << " Extens    
446 #  endif                                          
447                                                   
448   fMateIDs = new size_t[fNoVoxelsX * fNoVoxels    
449   for (G4int iz = 0; iz < fNoVoxelsZ; ++iz) {     
450     for (G4int iy = 0; iy < fNoVoxelsY; ++iy)     
451       for (G4int ix = 0; ix < fNoVoxelsX; ++ix    
452         G4int mateID;                             
453         fin >> mateID;                            
454         G4int nnew = ix + (iy)*fNoVoxelsX + (i    
455         if (mateID < 0 || mateID >= nMaterials    
456           G4Exception("GmReadPhantomG4Geometry    
457                       FatalException,             
458                       G4String("It should be b    
459                                + G4UIcommand::    
460                                + G4UIcommand::    
461                         .c_str());                
462         }                                         
463         fMateIDs[nnew] = mateID;                  
464       }                                           
465     }                                             
466   }                                               
467                                                   
468   ReadVoxelDensities(fin);                        
469                                                   
470   fin.close();                                    
471 #endif                                            
472 }                                                 
473                                                   
474 //....oooOO0OOooo........oooOO0OOooo........oo    
475 void DicomDetectorConstruction::ReadVoxelDensi    
476 {                                                 
477   G4String stemp;                                 
478   std::map<G4int, std::pair<G4double, G4double    
479   std::map<G4int, std::pair<G4double, G4double    
480   for (G4int ii = 0; ii < G4int(fPhantomMateri    
481     densiMinMax[ii] = std::pair<G4double, G4do    
482   }                                               
483                                                   
484   char* part = std::getenv("DICOM_CHANGE_MATER    
485   G4double densityDiff = -1.;                     
486   if (part) densityDiff = G4UIcommand::Convert    
487                                                   
488   std::map<G4int, G4double> densityDiffs;         
489   for (G4int ii = 0; ii < G4int(fPhantomMateri    
490     densityDiffs[ii] = densityDiff;  // curren    
491   }                                               
492   //  densityDiffs[0] = 0.0001; //air             
493                                                   
494   //--- Calculate the average material density    
495   std::map<std::pair<G4Material*, G4int>, matI    
496                                                   
497   //---- Read the material densities              
498   G4double dens;                                  
499   for (G4int iz = 0; iz < fNoVoxelsZ; ++iz) {     
500     for (G4int iy = 0; iy < fNoVoxelsY; ++iy)     
501       for (G4int ix = 0; ix < fNoVoxelsX; ++ix    
502         fin >> dens;                              
503         G4int copyNo = ix + (iy)*fNoVoxelsX +     
504                                                   
505         if (densityDiff != -1.) continue;         
506                                                   
507         //--- store the minimum and maximum de    
508         mpite = densiMinMax.find(G4int(fMateID    
509         if (dens < (*mpite).second.first) (*mp    
510         if (dens > (*mpite).second.second) (*m    
511         //--- Get material from original list     
512         G4int mateID = G4int(fMateIDs[copyNo])    
513         std::map<G4int, G4Material*>::const_it    
514                                                   
515         //--- Check if density is equal to the    
516         if (std::fabs(dens - (*imite).second->    
517           continue;                               
518                                                   
519         //--- Build material name with fPhanto    
520         G4int densityBin = (G4int(dens / densi    
521                                                   
522         G4String mateName = (*imite).second->G    
523         //--- Look if it is the first voxel wi    
524         std::pair<G4Material*, G4int> matdens(    
525                                                   
526         auto mppite = newMateDens.find(matdens    
527         if (mppite != newMateDens.cend()) {       
528           matInfo* mi = (*mppite).second;         
529           mi->fSumdens += dens;                   
530           mi->fNvoxels++;                         
531           fMateIDs[copyNo] = fPhantomMaterials    
532         }                                         
533         else {                                    
534           matInfo* mi = new matInfo;              
535           mi->fSumdens = dens;                    
536           mi->fNvoxels = 1;                       
537           mi->fId = G4int(newMateDens.size() +    
538           newMateDens[matdens] = mi;              
539           fMateIDs[copyNo] = fPhantomMaterials    
540         }                                         
541       }                                           
542     }                                             
543   }                                               
544                                                   
545   if (densityDiff != -1.) {                       
546     for (mpite = densiMinMax.begin(); mpite !=    
547 #ifdef G4VERBOSE                                  
548       G4cout << "DicomDetectorConstruction::Re    
549              << " ORIG MATERIALS DENSITY " <<     
550              << " MAX " << (*mpite).second.sec    
551 #endif                                            
552     }                                             
553   }                                               
554                                                   
555   //----- Build the list of phantom materials     
556   //--- Add original materials                    
557   for (auto mimite = fPhantomMaterialsOriginal    
558        ++mimite)                                  
559   {                                               
560     fMaterials.push_back((*mimite).second);       
561   }                                               
562   //                                              
563   //---- Build and add new materials              
564   for (auto mppite = newMateDens.cbegin(); mpp    
565     G4double averdens = (*mppite).second->fSum    
566     G4double saverdens = G4int(1000.001 * aver    
567 #ifndef G4VERBOSE                                 
568     G4cout << "DicomDetectorConstruction::Read    
569            << saverdens << " -> " << G4int(100    
570            << " " << G4int(1000 * averdens) /     
571 #endif                                            
572                                                   
573     G4String mateName =                           
574       ((*mppite).first).first->GetName() + "_"    
575     fMaterials.push_back(                         
576       BuildMaterialWithChangingDensity((*mppit    
577   }                                               
578 }                                                 
579                                                   
580 //....oooOO0OOooo........oooOO0OOooo........oo    
581 void DicomDetectorConstruction::ReadPhantomDat    
582 {                                                 
583   G4String dataFile = DicomHandler::GetDicomDa    
584   std::ifstream finDF(dataFile.c_str());          
585   G4String fname;                                 
586                                                   
587   if (finDF.good() != 1) {                        
588     G4String descript = "Problem reading data     
589     G4Exception(" DicomDetectorConstruction::R    
590   }                                               
591                                                   
592   G4int compression;                              
593   finDF >> compression;  // not used here         
594   finDF >> fNoFiles;                              
595                                                   
596   for (G4int i = 0; i < fNoFiles; ++i) {          
597     finDF >> fname;                               
598                                                   
599     //--- Read one data file                      
600     fname += ".g4dcm";                            
601                                                   
602     ReadPhantomDataFile(fname);                   
603   }                                               
604                                                   
605   //----- Merge data headers                      
606   MergeZSliceHeaders();                           
607   finDF.close();                                  
608 }                                                 
609                                                   
610 //....oooOO0OOooo........oooOO0OOooo........oo    
611 void DicomDetectorConstruction::ReadPhantomDat    
612 {                                                 
613   G4cout << " DicomDetectorConstruction::ReadP    
614          << G4endl;  // GDEB                      
615                                                   
616 #ifdef G4VERBOSE                                  
617   G4cout << " DicomDetectorConstruction::ReadP    
618 #endif                                            
619                                                   
620   std::ifstream fin(fname.c_str(), std::ios_ba    
621   if (!fin.is_open()) {                           
622     G4Exception("DicomDetectorConstruction::Re    
623                 G4String("File not found " + f    
624   }                                               
625   //----- Define density differences (maximum     
626   // a new material)                              
627   char* part = std::getenv("DICOM_CHANGE_MATER    
628   G4double densityDiff = -1.;                     
629   if (part) densityDiff = G4UIcommand::Convert    
630   if (densityDiff != -1.) {                       
631     for (unsigned int ii = 0; ii < fOriginalMa    
632       fDensityDiffs[ii] = densityDiff;  // cur    
633       // same difference                          
634     }                                             
635   }                                               
636   else {                                          
637     if (fMaterials.size() == 0) {  // do it on    
638       for (unsigned int ii = 0; ii < fOriginal    
639         fMaterials.push_back(fOriginalMaterial    
640       }                                           
641     }                                             
642   }                                               
643                                                   
644   //----- Read data header                        
645   DicomPhantomZSliceHeader* sliceHeader = new     
646   fZSliceHeaders.push_back(sliceHeader);          
647                                                   
648   //----- Read material indices                   
649   G4int nVoxels = sliceHeader->GetNoVoxels();     
650                                                   
651   //--- If first slice, initiliaze fMateIDs       
652   if (fZSliceHeaders.size() == 1) {               
653     // fMateIDs = new unsigned int[fNoFiles*nV    
654     fMateIDs = new size_t[fNoFiles * nVoxels];    
655   }                                               
656                                                   
657   unsigned int mateID;                            
658   // number of voxels from previously read sli    
659   G4int voxelCopyNo = G4int((fZSliceHeaders.si    
660   for (G4int ii = 0; ii < nVoxels; ++ii, voxel    
661     fin >> mateID;                                
662     fMateIDs[voxelCopyNo] = mateID;               
663   }                                               
664                                                   
665   //----- Read material densities and build ne    
666   //  same material but its density is in a di    
667   // (size of density intervals defined by den    
668   G4double density;                               
669   // number of voxels from previously read sli    
670   voxelCopyNo = G4int((fZSliceHeaders.size() -    
671   for (G4int ii = 0; ii < nVoxels; ++ii, voxel    
672     fin >> density;                               
673                                                   
674     //-- Get material from list of original ma    
675     mateID = unsigned(fMateIDs[voxelCopyNo]);     
676     G4Material* mateOrig = fOriginalMaterials[    
677                                                   
678     //-- Get density bin: middle point of the     
679     // density is included                        
680     G4String newMateName = mateOrig->GetName()    
681     G4float densityBin = 0.;                      
682     if (densityDiff != -1.) {                     
683       densityBin = G4float(fDensityDiffs[mateI    
684       //-- Build the new material name            
685       newMateName += G4UIcommand::ConvertToStr    
686     }                                             
687                                                   
688     //-- Look if a material with this name is     
689     //  (because a previous voxel was already     
690     unsigned int im;                              
691     for (im = 0; im < fMaterials.size(); ++im)    
692       if (fMaterials[im]->GetName() == newMate    
693         break;                                    
694       }                                           
695     }                                             
696     //-- If material is already created use in    
697     if (im != fMaterials.size()) {                
698       fMateIDs[voxelCopyNo] = im;                 
699       //-- else, create the material              
700     }                                             
701     else {                                        
702       if (densityDiff != -1.) {                   
703         fMaterials.push_back(BuildMaterialWith    
704         fMateIDs[voxelCopyNo] = fMaterials.siz    
705       }                                           
706       else {                                      
707         G4cerr << " im " << im << " < " << fMa    
708         G4Exception("DicomDetectorConstruction    
709                     "Wrong index in material")    
710       }                                           
711     }                                             
712   }                                               
713 }                                                 
714                                                   
715 //....oooOO0OOooo........oooOO0OOooo........oo    
716 void DicomDetectorConstruction::MergeZSliceHea    
717 {                                                 
718   //----- Images must have the same dimension     
719   fZSliceHeaderMerged = new DicomPhantomZSlice    
720   for (unsigned int ii = 1; ii < fZSliceHeader    
721     *fZSliceHeaderMerged += *fZSliceHeaders[ii    
722   }                                               
723 }                                                 
724                                                   
725 //....oooOO0OOooo........oooOO0OOooo........oo    
726 G4Material* DicomDetectorConstruction::BuildMa    
727                                                   
728                                                   
729 {                                                 
730   //----- Copy original material, but with new    
731   G4int nelem = G4int(origMate->GetNumberOfEle    
732   G4Material* mate =                              
733     new G4Material(newMateName, density * g /     
734                                                   
735   for (G4int ii = 0; ii < nelem; ++ii) {          
736     G4double frac = origMate->GetFractionVecto    
737     G4Element* elem = const_cast<G4Element*>(o    
738     mate->AddElement(elem, frac);                 
739   }                                               
740                                                   
741   return mate;                                    
742 }                                                 
743                                                   
744 //....oooOO0OOooo........oooOO0OOooo........oo    
745 void DicomDetectorConstruction::ConstructPhant    
746 {                                                 
747   //---- Extract number of voxels and voxel di    
748   fNoVoxelsX = fZSliceHeaderMerged->GetNoVoxel    
749   fNoVoxelsY = fZSliceHeaderMerged->GetNoVoxel    
750   fNoVoxelsZ = fZSliceHeaderMerged->GetNoVoxel    
751                                                   
752   fVoxelHalfDimX = fZSliceHeaderMerged->GetVox    
753   fVoxelHalfDimY = fZSliceHeaderMerged->GetVox    
754   fVoxelHalfDimZ = fZSliceHeaderMerged->GetVox    
755 #ifdef G4VERBOSE                                  
756   G4cout << " fNoVoxelsX " << fNoVoxelsX << "     
757   G4cout << " fNoVoxelsY " << fNoVoxelsY << "     
758   G4cout << " fNoVoxelsZ " << fNoVoxelsZ << "     
759   G4cout << " totalPixels " << fNoVoxelsX * fN    
760 #endif                                            
761                                                   
762   //----- Define the volume that contains all     
763   fContainer_solid = new G4Box("phantomContain    
764                                fNoVoxelsY * fV    
765   fContainer_logic =                              
766     new G4LogicalVolume(fContainer_solid,         
767                         // the material is not    
768                         fMaterials[0], "phanto    
769   //--- Place it on the world                     
770   G4double fOffsetX = (fZSliceHeaderMerged->Ge    
771   G4double fOffsetY = (fZSliceHeaderMerged->Ge    
772   G4double fOffsetZ = (fZSliceHeaderMerged->Ge    
773   G4ThreeVector posCentreVoxels(fOffsetX, fOff    
774 #ifdef G4VERBOSE                                  
775   G4cout << " placing voxel container volume a    
776 #endif                                            
777   fContainer_phys = new G4PVPlacement(0,  // r    
778                                       posCentr    
779                                       fContain    
780                                       "phantom    
781                                       fWorld_l    
782                                       false,      
783                                       1);  //     
784 }                                                 
785                                                   
786 //....oooOO0OOooo........oooOO0OOooo........oo    
787 void DicomDetectorConstruction::ConstructPhant    
788 {                                                 
789 #ifdef G4_DCMTK                                   
790   //---- Extract number of voxels and voxel di    
791 #  ifdef G4VERBOSE                                
792   G4cout << " fNoVoxelsX " << fNoVoxelsX << "     
793   G4cout << " fNoVoxelsY " << fNoVoxelsY << "     
794   G4cout << " fNoVoxelsZ " << fNoVoxelsZ << "     
795   G4cout << " totalPixels " << fNoVoxelsX * fN    
796 #  endif                                          
797                                                   
798   //----- Define the volume that contains all     
799   fContainer_solid = new G4Box("phantomContain    
800                                fNoVoxelsY * fV    
801   fContainer_logic =                              
802     new G4LogicalVolume(fContainer_solid,         
803                         // the material is not    
804                         fMaterials[0], "phanto    
805                                                   
806   G4ThreeVector posCentreVoxels((fMinX + fMaxX    
807 #  ifdef G4VERBOSE                                
808   G4cout << " placing voxel container volume a    
809 #  endif                                          
810   fContainer_phys = new G4PVPlacement(0,  // r    
811                                       posCentr    
812                                       fContain    
813                                       "phantom    
814                                       fWorld_l    
815                                       false,      
816                                       1);  //     
817 #endif                                            
818 }                                                 
819                                                   
820 #include "G4MultiFunctionalDetector.hh"           
821 #include "G4PSDoseDeposit.hh"                     
822 #include "G4PSDoseDeposit3D.hh"                   
823 #include "G4SDManager.hh"                         
824                                                   
825 //....oooOO0OOooo........oooOO0OOooo........oo    
826 void DicomDetectorConstruction::SetScorer(G4Lo    
827 {                                                 
828 #ifdef G4VERBOSE                                  
829   G4cout << "\t SET SCORER : " << voxel_logic-    
830 #endif                                            
831                                                   
832   fScorers.insert(voxel_logic);                   
833 }                                                 
834                                                   
835 //....oooOO0OOooo........oooOO0OOooo........oo    
836                                                   
837 void DicomDetectorConstruction::ConstructSDand    
838 {                                                 
839 #ifdef G4VERBOSE                                  
840   G4cout << "\t CONSTRUCT SD AND FIELD" << G4e    
841 #endif                                            
842                                                   
843   // G4SDManager* SDman = G4SDManager::GetSDMp    
844                                                   
845   // SDman->SetVerboseLevel(1);                   
846                                                   
847   //                                              
848   // Sensitive Detector Name                      
849   G4String concreteSDname = "phantomSD";          
850   std::vector<G4String> scorer_names;             
851   scorer_names.push_back(concreteSDname);         
852   //------------------------                      
853   // MultiFunctionalDetector                      
854   //------------------------                      
855   //                                              
856   // Define MultiFunctionalDetector with name.    
857   // declare MFDet as a MultiFunctionalDetecto    
858   G4MultiFunctionalDetector* MFDet = new G4Mul    
859   G4SDManager::GetSDMpointer()->AddNewDetector    
860   char* nest = std::getenv("DICOM_NESTED_PARAM    
861   if (nest && G4String(nest) == "1") {            
862     G4VPrimitiveScorer* dosedep = new G4PSDose    
863       "DoseDeposit", fNoVoxelsZ, fNoVoxelsY, f    
864     // - the last 3 arguments correspond to th    
865     MFDet->RegisterPrimitive(dosedep);            
866   }                                               
867   else {                                          
868     G4VPrimitiveScorer* dosedep = new G4PSDose    
869     MFDet->RegisterPrimitive(dosedep);            
870   }                                               
871                                                   
872   for (auto ite = fScorers.cbegin(); ite != fS    
873     SetSensitiveDetector(*ite, MFDet);            
874   }                                               
875 }                                                 
876