Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/ICRP110_HumanPhantoms/src/ICRP110PhantomConstruction.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 ]

  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 // Code developed by:
 27 // S.Guatelli, M. Large and A. Malaroda, University of Wollongong
 28 //
 29 #include "ICRP110PhantomConstruction.hh"
 30 #include "ICRP110PhantomNestedParameterisation.hh"
 31 #include "ICRP110PhantomMaterial_Female.hh"
 32 #include "ICRP110PhantomMaterial_Male.hh"
 33 #include "globals.hh"
 34 #include "G4SystemOfUnits.hh"
 35 #include "G4Box.hh"
 36 #include "G4Colour.hh"
 37 #include "G4LogicalVolume.hh"
 38 #include "G4PVPlacement.hh"
 39 #include "G4PVParameterised.hh"
 40 #include "G4VPhysicalVolume.hh"
 41 #include "G4PVPlacement.hh"
 42 #include "G4PVParameterised.hh"
 43 #include "G4RunManager.hh"
 44 #include "G4VisAttributes.hh"
 45 #include <map>
 46 #include <cstdlib>
 47 
 48 ICRP110PhantomConstruction::ICRP110PhantomConstruction():
 49    fMotherVolume(nullptr), fPhantomContainer(nullptr),
 50    fNVoxelX(0), fNVoxelY(0), fNVoxelZ(0), 
 51    fVoxelHalfDimX(0), fVoxelHalfDimY(0), fVoxelHalfDimZ(0),
 52    fMinX(0),fMaxX(0), fMinY(0), fMaxY(0),
 53    fMinZ(0), fMaxZ(0), fNoFiles(0), fNVoxels(0),
 54    fMateIDs(nullptr)
 55 {
 56   fMessenger = new ICRP110PhantomMessenger(this);
 57   // the messenger allows to set the sex of the phantom
 58   // interactively
 59   fMaterial_Female = new ICRP110PhantomMaterial_Female();
 60   fMaterial_Male = new ICRP110PhantomMaterial_Male();
 61   fSex = "female"; // Female phantom is the default option
 62   fSection = "head"; // Head partial phantom is the default option
 63 }
 64 
 65 ICRP110PhantomConstruction::~ICRP110PhantomConstruction()
 66 {
 67   delete fMaterial_Female;
 68   delete fMaterial_Male;
 69   delete fMessenger;
 70 }
 71 
 72 G4VPhysicalVolume* ICRP110PhantomConstruction::Construct()
 73 {
 74   // Define Material Air
 75   G4double A;  // atomic mass
 76   G4double Z;  // atomic number
 77   G4double d;  // density
 78   
 79   A = 14.01*g/mole;
 80   auto elN = new G4Element("Nitrogen","N",Z = 7.,A);
 81   A = 16.00*g/mole;
 82   auto elO = new G4Element("Oxygen","O",Z = 8.,A);
 83   
 84   d = 0.001 *g/cm3;
 85   auto matAir = new G4Material("Air",d,2);
 86   matAir -> AddElement(elN,0.8);
 87   matAir -> AddElement(elO,0.2); 
 88 
 89  std::vector<G4Material*> pMaterials;
 90 
 91  if(fSex == "female"){
 92 
 93   fMaterial_Female -> DefineMaterials();
 94 //----- Store materials in a vector
 95     pMaterials.push_back(matAir); 
 96     pMaterials.push_back(fMaterial_Female -> GetMaterial("teeth")); 
 97     pMaterials.push_back(fMaterial_Female -> GetMaterial("bone"));
 98     pMaterials.push_back(fMaterial_Female -> GetMaterial("humeri_upper")); 
 99     pMaterials.push_back(fMaterial_Female -> GetMaterial("humeri_lower")); 
100     pMaterials.push_back(fMaterial_Female -> GetMaterial("arm_lower")); 
101     pMaterials.push_back(fMaterial_Female -> GetMaterial("hand")); 
102     pMaterials.push_back(fMaterial_Female -> GetMaterial("clavicle")); 
103     pMaterials.push_back(fMaterial_Female -> GetMaterial("cranium")); 
104     pMaterials.push_back(fMaterial_Female -> GetMaterial("femora_upper")); 
105     pMaterials.push_back(fMaterial_Female -> GetMaterial("femora_lower")); 
106     pMaterials.push_back(fMaterial_Female -> GetMaterial("leg_lower")); 
107     pMaterials.push_back(fMaterial_Female -> GetMaterial("foot")); 
108     pMaterials.push_back(fMaterial_Female -> GetMaterial("mandible")); 
109     pMaterials.push_back(fMaterial_Female -> GetMaterial("pelvis")); 
110     pMaterials.push_back(fMaterial_Female -> GetMaterial("ribs")); 
111     pMaterials.push_back(fMaterial_Female -> GetMaterial("scapulae")); 
112     pMaterials.push_back(fMaterial_Female -> GetMaterial("spine_cervical")); 
113     pMaterials.push_back(fMaterial_Female -> GetMaterial("spine_thoratic")); 
114     pMaterials.push_back(fMaterial_Female -> GetMaterial("spine_lumbar")); 
115     pMaterials.push_back(fMaterial_Female -> GetMaterial("sacrum")); 
116     pMaterials.push_back(fMaterial_Female -> GetMaterial("sternum")); 
117     pMaterials.push_back(fMaterial_Female -> GetMaterial("hf_upper")); 
118     pMaterials.push_back(fMaterial_Female -> GetMaterial("hf_lower")); 
119     pMaterials.push_back(fMaterial_Female -> GetMaterial("med_lowerarm")); 
120     pMaterials.push_back(fMaterial_Female -> GetMaterial("med_lowerleg")); 
121     pMaterials.push_back(fMaterial_Female -> GetMaterial("cartilage")); 
122     pMaterials.push_back(fMaterial_Female -> GetMaterial("skin")); 
123     pMaterials.push_back(fMaterial_Female -> GetMaterial("blood")); 
124     pMaterials.push_back(fMaterial_Female -> GetMaterial("muscle")); 
125     pMaterials.push_back(fMaterial_Female -> GetMaterial("liver")); 
126     pMaterials.push_back(fMaterial_Female -> GetMaterial("pancreas")); 
127     pMaterials.push_back(fMaterial_Female -> GetMaterial("brain"));
128     pMaterials.push_back(fMaterial_Female -> GetMaterial("heart")); 
129     pMaterials.push_back(fMaterial_Female -> GetMaterial("eye")); 
130     pMaterials.push_back(fMaterial_Female -> GetMaterial("kidney")); 
131     pMaterials.push_back(fMaterial_Female -> GetMaterial("stomach")); 
132     pMaterials.push_back(fMaterial_Female -> GetMaterial("intestine_sml")); 
133     pMaterials.push_back(fMaterial_Female -> GetMaterial("intestine_lrg")); 
134     pMaterials.push_back(fMaterial_Female -> GetMaterial("spleen")); 
135     pMaterials.push_back(fMaterial_Female -> GetMaterial("thyroid")); 
136     pMaterials.push_back(fMaterial_Female -> GetMaterial("bladder")); 
137     pMaterials.push_back(fMaterial_Female -> GetMaterial("ovaries_testes")); 
138     pMaterials.push_back(fMaterial_Female -> GetMaterial("adrenals")); 
139     pMaterials.push_back(fMaterial_Female -> GetMaterial("oesophagus")); 
140     pMaterials.push_back(fMaterial_Female -> GetMaterial("misc")); 
141     pMaterials.push_back(fMaterial_Female -> GetMaterial("uterus_prostate"));
142     pMaterials.push_back(fMaterial_Female -> GetMaterial("lymph")); 
143     pMaterials.push_back(fMaterial_Female -> GetMaterial("breast_glandular")); 
144     pMaterials.push_back(fMaterial_Female -> GetMaterial("breast_adipose")); 
145     pMaterials.push_back(fMaterial_Female -> GetMaterial("lung")); 
146     pMaterials.push_back(fMaterial_Female -> GetMaterial("gastro_content")); 
147     pMaterials.push_back(fMaterial_Female -> GetMaterial("urine"));     
148  }
149  else if (fSex == "male"){
150  // MATT do the same here
151     fMaterial_Male -> DefineMaterials();
152 
153 //----- Store materials in a vector
154     pMaterials.push_back(matAir); 
155     pMaterials.push_back(fMaterial_Male -> GetMaterial("teeth")); 
156     pMaterials.push_back(fMaterial_Male -> GetMaterial("bone"));
157     pMaterials.push_back(fMaterial_Male -> GetMaterial("humeri_upper")); 
158     pMaterials.push_back(fMaterial_Male -> GetMaterial("humeri_lower")); 
159     pMaterials.push_back(fMaterial_Male -> GetMaterial("arm_lower")); 
160     pMaterials.push_back(fMaterial_Male -> GetMaterial("hand")); 
161     pMaterials.push_back(fMaterial_Male -> GetMaterial("clavicle")); 
162     pMaterials.push_back(fMaterial_Male -> GetMaterial("cranium")); 
163     pMaterials.push_back(fMaterial_Male -> GetMaterial("femora_upper")); 
164     pMaterials.push_back(fMaterial_Male -> GetMaterial("femora_lower")); 
165     pMaterials.push_back(fMaterial_Male -> GetMaterial("leg_lower")); 
166     pMaterials.push_back(fMaterial_Male -> GetMaterial("foot")); 
167     pMaterials.push_back(fMaterial_Male -> GetMaterial("mandible")); 
168     pMaterials.push_back(fMaterial_Male -> GetMaterial("pelvis")); 
169     pMaterials.push_back(fMaterial_Male -> GetMaterial("ribs")); 
170     pMaterials.push_back(fMaterial_Male -> GetMaterial("scapulae")); 
171     pMaterials.push_back(fMaterial_Male -> GetMaterial("spine_cervical")); 
172     pMaterials.push_back(fMaterial_Male -> GetMaterial("spine_thoratic")); 
173     pMaterials.push_back(fMaterial_Male -> GetMaterial("spine_lumbar")); 
174     pMaterials.push_back(fMaterial_Male -> GetMaterial("sacrum")); 
175     pMaterials.push_back(fMaterial_Male -> GetMaterial("sternum")); 
176     pMaterials.push_back(fMaterial_Male -> GetMaterial("hf_upper")); 
177     pMaterials.push_back(fMaterial_Male -> GetMaterial("hf_lower")); 
178     pMaterials.push_back(fMaterial_Male -> GetMaterial("med_lowerarm")); 
179     pMaterials.push_back(fMaterial_Male -> GetMaterial("med_lowerleg")); 
180     pMaterials.push_back(fMaterial_Male -> GetMaterial("cartilage")); 
181     pMaterials.push_back(fMaterial_Male -> GetMaterial("skin")); 
182     pMaterials.push_back(fMaterial_Male -> GetMaterial("blood")); 
183     pMaterials.push_back(fMaterial_Male -> GetMaterial("muscle")); 
184     pMaterials.push_back(fMaterial_Male -> GetMaterial("liver")); 
185     pMaterials.push_back(fMaterial_Male -> GetMaterial("pancreas")); 
186     pMaterials.push_back(fMaterial_Male -> GetMaterial("brain"));
187     pMaterials.push_back(fMaterial_Male -> GetMaterial("heart")); 
188     pMaterials.push_back(fMaterial_Male -> GetMaterial("eye")); 
189     pMaterials.push_back(fMaterial_Male -> GetMaterial("kidney")); 
190     pMaterials.push_back(fMaterial_Male -> GetMaterial("stomach")); 
191     pMaterials.push_back(fMaterial_Male -> GetMaterial("intestine_sml")); 
192     pMaterials.push_back(fMaterial_Male -> GetMaterial("intestine_lrg")); 
193     pMaterials.push_back(fMaterial_Male -> GetMaterial("spleen")); 
194     pMaterials.push_back(fMaterial_Male -> GetMaterial("thyroid")); 
195     pMaterials.push_back(fMaterial_Male -> GetMaterial("bladder")); 
196     pMaterials.push_back(fMaterial_Male -> GetMaterial("ovaries_testes")); 
197     pMaterials.push_back(fMaterial_Male -> GetMaterial("adrenals")); 
198     pMaterials.push_back(fMaterial_Male -> GetMaterial("oesophagus")); 
199     pMaterials.push_back(fMaterial_Male -> GetMaterial("misc")); 
200     pMaterials.push_back(fMaterial_Male -> GetMaterial("uterus_prostate"));
201     pMaterials.push_back(fMaterial_Male -> GetMaterial("lymph")); 
202     pMaterials.push_back(fMaterial_Male -> GetMaterial("breast_glandular")); 
203     pMaterials.push_back(fMaterial_Male -> GetMaterial("breast_adipose")); 
204     pMaterials.push_back(fMaterial_Male -> GetMaterial("lung")); 
205     pMaterials.push_back(fMaterial_Male -> GetMaterial("gastro_content")); 
206     pMaterials.push_back(fMaterial_Male -> GetMaterial("urine")); 
207  
208 }
209   
210   // World Volume
211   G4double worldSize = 2.*m ;
212   G4Box* world = new G4Box("world", worldSize, worldSize, worldSize);
213 
214   auto logicWorld = new G4LogicalVolume(world,
215                        matAir,
216                  "logicalWorld", nullptr, nullptr,nullptr);
217 
218   fMotherVolume = new G4PVPlacement(nullptr,G4ThreeVector(),
219             "physicalWorld",
220             logicWorld,
221             nullptr,
222             false,
223             0);
224 
225   logicWorld -> SetVisAttributes(G4VisAttributes::GetInvisible());
226  
227   G4cout << "World has been built" << G4endl; 
228 
229   G4cout << "Phantom Sex: " << fSex << G4endl;
230   G4cout << "Phantom Section: " << fSection << G4endl;
231   ReadPhantomData(fSex, fSection); 
232   
233   G4cout << "Number of X,Y,Z voxels = " << fNVoxelX << ", " << fNVoxelY << ", " << fNVoxelZ << G4endl;
234   
235 //----- Define the volume that contains all the voxels
236   G4Box* fContainer_solid = new G4Box("phantomContainer",fNVoxelX*fVoxelHalfDimX*mm,
237                                fNVoxelY*fVoxelHalfDimY*mm,
238                                fNVoxelZ*fVoxelHalfDimZ*mm);
239  
240   auto fContainer_logic = new G4LogicalVolume( fContainer_solid,
241                                                             matAir,
242                                                             "phantomContainer",
243                                                              nullptr, nullptr, nullptr);                                                        
244     
245   fMaxX = fNVoxelX*fVoxelHalfDimX*mm; // Max X along X axis of the voxelised geometry 
246   fMaxY = fNVoxelY*fVoxelHalfDimY*mm; // Max Y
247   fMaxZ = fNVoxelZ*fVoxelHalfDimZ*mm; // Max Z
248 
249   fMinX = -fNVoxelX*fVoxelHalfDimX*mm;// Min X 
250   fMinY = -fNVoxelY*fVoxelHalfDimY*mm;// Min Y
251   fMinZ = -fNVoxelZ*fVoxelHalfDimZ*mm;// Min Z
252 
253   G4ThreeVector posCentreVoxels((fMinX+fMaxX)/2.,(fMinY+fMaxY)/2.,(fMinZ+fMaxZ)/2.);
254 
255   G4cout << " placing voxel container volume at " << posCentreVoxels << G4endl;
256 
257    
258   fPhantomContainer
259   = new G4PVPlacement(nullptr,                     // rotation
260                       posCentreVoxels,
261                       fContainer_logic,     // The logic volume
262                       "phantomContainer",  // Name
263                       logicWorld,         // Mother
264                       false,            // No op. bool.
265                       1);              // Copy number
266   
267   fContainer_logic -> SetVisAttributes(new G4VisAttributes(G4Colour(1.,0.,0.,0.)));
268 
269 
270 // Define the voxelised phantom here
271 // Replication of air Phantom Volume.
272 
273 //--- Slice the phantom along Y axis
274    G4String yRepName("RepY");
275    G4VSolid* solYRep = new G4Box(yRepName,fNVoxelX*fVoxelHalfDimX,
276                                   fVoxelHalfDimY, fNVoxelZ*fVoxelHalfDimZ);
277    auto logYRep = new G4LogicalVolume(solYRep,matAir,yRepName);
278    new G4PVReplica(yRepName,logYRep,fContainer_logic,kYAxis, fNVoxelY,fVoxelHalfDimY*2.);
279   
280    logYRep -> SetVisAttributes(new G4VisAttributes(G4VisAttributes::GetInvisible()));   
281 
282 //--- Slice the phantom along X axis 
283    G4String xRepName("RepX");
284    G4VSolid* solXRep = new G4Box(xRepName,fVoxelHalfDimX,fVoxelHalfDimY,
285                                   fNVoxelZ*fVoxelHalfDimZ);
286    auto logXRep = new G4LogicalVolume(solXRep,matAir,xRepName);
287    new G4PVReplica(xRepName,logXRep,logYRep,kXAxis,fNVoxelX,fVoxelHalfDimX*2.);
288 
289    logXRep -> SetVisAttributes(new G4VisAttributes(G4VisAttributes::GetInvisible()));
290   
291    //----- Voxel solid and logical volumes
292    //--- Slice along Z axis 
293    G4VSolid* solidVoxel = new G4Box("phantom",fVoxelHalfDimX, fVoxelHalfDimY,fVoxelHalfDimZ);
294    auto logicVoxel = new G4LogicalVolume(solidVoxel,matAir,"phantom");
295 
296    logicVoxel -> SetVisAttributes(new G4VisAttributes(G4VisAttributes::GetInvisible()));
297 
298    // Parameterisation to define the material of each voxel
299    G4ThreeVector halfVoxelSize(fVoxelHalfDimX,fVoxelHalfDimY,fVoxelHalfDimZ);
300       
301    auto param =  new ICRP110PhantomNestedParameterisation(halfVoxelSize, pMaterials);
302 
303    new G4PVParameterised("phantom",    // their name
304                           logicVoxel, // their logical volume
305                           logXRep,      // Mother logical volume
306                           kZAxis,       // Are placed along this axis
307                           fNVoxelZ,      // Number of cells
308                           param);       // Parameterisation
309 
310     param -> SetMaterialIndices(fMateIDs); // fMateIDs is  the vector with Material ID associated to each voxel, from ASCII input data files.
311     param -> SetNoVoxel(fNVoxelX,fNVoxelY,fNVoxelZ);
312 
313   return fMotherVolume;
314 }
315 
316 void ICRP110PhantomConstruction::ReadPhantomData(const G4String& sex, const G4String& section)
317 {
318 
319   // This method reads the information of ICRPdata/FemaleData.dat or
320   // ICRPdata/MaleData.data depending on the sex of the chosen phantom
321 
322 fSex = sex;
323 fSection = section;
324 
325 G4String dataFile;
326 
327     if (fSex == "female")
328     {
329         if (fSection == "head")
330         {
331           dataFile = "ICRPdata/FemaleHead.dat";
332         }
333         else if (fSection == "trunk")
334         {
335           dataFile = "ICRPdata/FemaleTrunk.dat";
336         }
337         else if (fSection == "full")
338         {
339           dataFile = "ICRPdata/FemaleData.dat"; 
340         }
341     }
342     if (fSex == "male")
343     {
344         if (fSection == "head")
345         {
346           dataFile = "ICRPdata/MaleHead.dat";
347         }
348         else if (fSection == "trunk")
349         {
350           dataFile = "ICRPdata/MaleTrunk.dat";
351         }
352         else if (fSection == "full")
353         {
354           dataFile = "ICRPdata/MaleData.dat"; 
355         }
356     }
357     
358     G4cout << "Data file " << dataFile << " is read by Detector Construction." << G4endl;
359     
360 // The data.dat file in directory/build/ICRPdata/ contains the information 
361 // to build the phantoms. For more details look in the README file.
362 
363   //input file named finDF which consists of dataFile as a string object
364   std::ifstream finDF(dataFile.c_str()); 
365   
366  
367   G4String fname;
368 
369 if(finDF.good() != 1 ) //check that the file is good and working
370  { 
371   G4String descript = "Problem reading data file: "+dataFile;
372   G4Exception(" HumanPhantomConstruction::ReadPhantomData"," ", 
373               FatalException,descript);
374   }
375 
376   finDF >> fNoFiles;
377   G4cout << "Number of files = " << fNoFiles << G4endl;
378   finDF >> fNVoxelX;      //Inputs number of X-Voxels
379   finDF >> fNVoxelY;      //Y-Voxels
380   fNVoxelZ = fNoFiles;    //Z-Voxels (equal to number of slice files built/read)
381   finDF >> fVoxelHalfDimX;
382   finDF >> fVoxelHalfDimY;
383   finDF >> fVoxelHalfDimZ;
384   G4cout << "Number of X,Y,Z voxels = " << fNVoxelX << ", " << fNVoxelY << ", " << fNVoxelZ <<G4endl;
385   
386   fNVoxels = fNVoxelX*fNVoxelY*fNVoxelZ; 
387   G4cout << "Total Number of Voxels = " << fNVoxels << G4endl;
388   
389   G4int nMaterials;
390   finDF >> nMaterials;
391   G4String mateName;
392   G4int nmate;
393 
394 //-----Read materials and associate with material ID number------//
395   
396   for( G4int ii = 0; ii < nMaterials; ii++ ){
397     finDF >> nmate;
398     finDF >> mateName;
399     
400     // This allows to skip empty spaces and tabs in the string 
401     if( mateName[0] == '"' && mateName[mateName.length()-1] == '"' ) 
402     {
403       mateName = mateName.substr(1,mateName.length()-2); 
404     }
405  
406     // To uncomment for eventual debugging
407     /* G4cout << "GmReadPhantomG4Geometry::ReadPhantomData reading nmate " 
408            << ii << " = " << nmate 
409            << " mate " << mateName << G4endl;*/
410  
411     if( ii != nmate ) {
412     G4Exception("GmReadPhantomG4Geometry::ReadPhantomData",
413                 "Wrong argument",
414                 FatalErrorInArgument,
415                 "Material number should be in increasing order:wrong material number");
416                 }
417       }
418   
419 fMateIDs = new size_t[fNVoxels]; //Array with Material ID for each voxel
420 
421 G4cout << "ICRP110PhantomConstruction::ReadPhantomDataFile is openining the following phantom files: " << G4endl;
422   
423 for(G4int i = 0; i < fNoFiles; i++ )
424   {
425     finDF >> fname;
426     ReadPhantomDataFile(fSex, fname, i); 
427   }
428 
429 finDF.close();
430 }
431 
432 //----------------Opens phantom ASCII slice files to construct the phantom from-----------------//
433   
434 void ICRP110PhantomConstruction::ReadPhantomDataFile(const G4String& sex, const G4String& fileName, G4int numberFile)
435 {
436 G4cout << fileName << G4endl;
437          
438 fSex = sex;
439 
440 G4String slice;
441 
442     if (fSex == "female")
443     {
444       slice = "ICRPdata/ICRP110_g4dat/AF/"+fileName;
445     }
446     if (fSex == "male")
447     {
448       slice = "ICRPdata/ICRP110_g4dat/AM/"+fileName;
449     }  
450   
451   std::ifstream fin(slice.c_str(), std::ios_base::in);
452   
453   if( !fin.is_open() ) {
454     G4Exception("HumanPhantomConstruction::ReadPhantomDataFile",
455                 "",
456                 FatalErrorInArgument,
457                 G4String("File not found " + fileName ).c_str());
458   }
459 
460     for( G4int iy = 0; iy < fNVoxelY; iy++ ) {
461       for( G4int ix = 0; ix < fNVoxelX; ix++ ) {
462       if (ix == 0 && iy == 0)
463         {
464           G4int dudX,dudY,dudZ;      
465           fin >> dudX >> dudY >> dudZ ;
466           // Dummy method to skip the first three lines of the files
467           // which are not used here
468         }
469      
470         G4int nnew = ix + (iy)*fNVoxelX + numberFile*fNVoxelX*fNVoxelY;
471         G4int OrgID;
472         fin >> OrgID; 
473 
474         G4int mateID_out;
475 
476 // The code below associates organ ID numbers (called here mateID) from ASCII slice
477 // files with material ID numbers (called here mateID_out) as defined in ICRP110PhantomMaterials
478 // Material and Organ IDs are associated as stated in AM_organs.dat and FM_organs.dat depending on
479 // the sex of the phantom (male and female, respctively)
480 
481   if (OrgID==128)
482   {
483   mateID_out=1;
484   }
485   
486   else if (OrgID==13 || OrgID==16 || OrgID==19 || OrgID==22 || OrgID==24 || OrgID==26 || OrgID==28 || OrgID==31 || OrgID==34 || OrgID==37 || OrgID==39 || OrgID==41 || OrgID==43 || OrgID==45 || OrgID==47 || OrgID==49 || OrgID==51 || OrgID==53 || OrgID==55)
487   {
488   mateID_out=2;
489   }
490 
491   else if (OrgID==14)
492   {
493   mateID_out=3;
494   }
495 
496   else if (OrgID==17)
497   {
498   mateID_out=4;
499   }
500 
501   else if (OrgID==20)
502   {
503   mateID_out=5;
504   }
505 
506   else if (OrgID==23)
507   {
508   mateID_out=6;
509   }
510 
511   else if (OrgID==25)
512   {
513   mateID_out=7;
514   }
515 
516   else if (OrgID==27)
517   {
518   mateID_out=8;
519   }
520 
521   else if (OrgID==29)
522   {
523   mateID_out=9;
524   }
525 
526   else if (OrgID==32)
527   {
528   mateID_out=10;
529   }
530 
531   else if (OrgID==35)
532   {
533   mateID_out=11;
534   }
535    
536   else if (OrgID==38)
537   {
538   mateID_out=12;
539   }
540 
541   else if (OrgID==40)
542   {
543   mateID_out=13;
544   }
545 
546   else if (OrgID==42)
547   {
548   mateID_out=14;
549   }
550 
551   else if (OrgID==44)
552   { 
553   mateID_out=15;
554   }
555 
556   else if (OrgID==46)
557   {
558   mateID_out=16;
559   }
560 
561   else if (OrgID==48)
562   {
563   mateID_out=17;
564   }
565 
566   else if (OrgID==50)
567   {
568   mateID_out=18;
569   }
570 
571   else if (OrgID==52)
572   {
573   mateID_out=19;
574   }
575 
576   else if (OrgID==54)
577   {
578   mateID_out=20;
579   }
580 
581   else if (OrgID==56)
582   {
583   mateID_out=21;
584   }
585 
586   else if (OrgID==15 || OrgID==30)
587   {
588   mateID_out=22;
589   }
590 
591   else if (OrgID==18 || OrgID==33)
592   {
593   mateID_out=23;
594   }
595 
596   else if (OrgID==21)
597   {
598   mateID_out=24;
599   }
600 
601   else if (OrgID==36)
602   {
603   mateID_out=25;
604   }
605 
606   else if (OrgID==57 || OrgID==58 || OrgID==59 || OrgID==60)  
607   {
608   mateID_out=26;
609   }
610 
611   else if (OrgID==122 || OrgID==123 || OrgID==124 || OrgID==125 || OrgID==141 )   
612   {
613   mateID_out=27;
614   }
615 
616   else if (OrgID==9 || OrgID==10 || OrgID==11 || OrgID==12 || OrgID==88 || OrgID==96 || OrgID==98)
617   {
618   mateID_out=28;
619   }
620 
621   else if (OrgID==5 || OrgID==6 || OrgID==106 || OrgID==107 || OrgID==108 || OrgID==109 || OrgID==133)
622   {
623   mateID_out=29;
624   }
625 
626   else if (OrgID==95)
627   {
628   mateID_out=30;
629   }
630 
631   else if (OrgID==113)
632   {
633   mateID_out=31;
634   }
635 
636   else if (OrgID==61)
637   {
638   mateID_out=32;
639   }
640 
641   else if (OrgID==87)
642   {
643   mateID_out=33;
644   }
645 
646   else if (OrgID==66 || OrgID==67 || OrgID==68 || OrgID==69)
647   {
648   mateID_out=34;
649   }
650 
651   else if (OrgID==89 || OrgID==90 || OrgID==91 || OrgID==92 || OrgID==93 || OrgID==94)
652   {
653   mateID_out=35;
654   }
655 
656   else if (OrgID==72)
657   {
658   mateID_out=36;
659   }
660 
661   else if (OrgID==74)
662   {
663   mateID_out=37;
664   }
665 
666   else if (OrgID==76 || OrgID==78 || OrgID==80 || OrgID==82 || OrgID==84 || OrgID==86)
667   {
668   mateID_out=38;
669   }
670 
671   else if (OrgID==127)
672   {
673   mateID_out=39;
674   }
675 
676   else if (OrgID==132)
677   {
678   mateID_out=40;
679   }
680 
681   else if (OrgID==137)
682   {
683   mateID_out=41;
684   }
685 
686   else if (OrgID==111 || OrgID==112 || OrgID==129 || OrgID==130)
687   {
688   mateID_out=42;
689   }
690 
691   else if (OrgID==1 || OrgID==2)
692   {
693   mateID_out=43;
694   }
695 
696   else if (OrgID==110)
697   {
698   mateID_out=44;
699   }
700 
701   else if (OrgID==3 || OrgID==4 || OrgID==7 || OrgID==8 || OrgID==70 || OrgID==71 || OrgID==114 || OrgID==120 || OrgID==121 || OrgID==126 || OrgID==131 || OrgID==134 || OrgID==135 || OrgID == 136)
702   {
703   mateID_out=45;
704   }
705 
706   else if (OrgID==115 || OrgID==139)
707   {
708   mateID_out=46;
709   }
710 
711   else if (OrgID==100 || OrgID==101 || OrgID==102 || OrgID==103 || OrgID==104 || OrgID==105)
712   {
713   mateID_out=47;
714   }
715 
716   else if (OrgID==63 || OrgID==65)
717   {
718   mateID_out=48;
719   }
720 
721   else if (OrgID==62 || OrgID==64 || OrgID==116 || OrgID==117 || OrgID==118 || OrgID==119)
722   {
723   mateID_out=49;
724   }
725 
726   else if (OrgID==97 || OrgID==99)
727   {
728   mateID_out=50;
729   }
730 
731   else if (OrgID==73 || OrgID==75 || OrgID==77 || OrgID==79 || OrgID==81 || OrgID==83 || OrgID==85)
732   {
733   mateID_out=51;
734   }
735 
736   else if (OrgID==138)
737   {
738   mateID_out=52;
739   }
740 
741   else if (OrgID==0 || OrgID==140)
742   {
743   mateID_out=0;
744   }
745 
746   else 
747   {
748   mateID_out=OrgID;
749   }
750         
751         G4int nMaterials = 53;
752         if( mateID_out < 0 || mateID_out >= nMaterials ) {
753           G4Exception("GmReadPhantomG4Geometry::ReadPhantomData",
754                       "Wrong index in phantom file",
755                       FatalException,
756                       G4String("It should be between 0 and "
757                               + G4UIcommand::ConvertToString(nMaterials-1) 
758                               + ", while it is " 
759                               + G4UIcommand::ConvertToString(OrgID)).c_str());
760         
761 //-------------Store Material IDs and position/reference number within phantom in vector---------------//
762    }
763   
764           fMateIDs[nnew] = mateID_out;
765   
766          
767       }
768    }
769 }
770 
771 //-----------Define phantom sex (male or female) to be constructed-------------//
772 void ICRP110PhantomConstruction::SetPhantomSex(G4String newSex)
773 {
774   fSex = newSex;
775   
776   if (fSex == "male")
777     {
778       G4cout << ">> Male Phantom will be built." << G4endl;
779     }
780   if (fSex == "female")
781     {
782       G4cout << ">> Female Phantom will be built." << G4endl;
783     }
784   if ((fSex != "female") && (fSex != "male"))
785     G4cout << fSex << " is not defined!" << G4endl;
786 } 
787   
788 //-----------Define phantom section to be constructed-------------//
789 void ICRP110PhantomConstruction::SetPhantomSection(G4String newSection)
790 {
791   fSection = newSection;
792   if (fSection == "head")
793     {
794       G4cout << ">> Partial Head Phantom will be built." << G4endl;
795     }
796   if (fSection == "trunk")
797     {
798       G4cout << ">> Partial Trunk Phantom will be built." << G4endl;
799     }
800   if (fSection == "full")
801     {
802       G4cout << ">> Full/Custom Phantom will be built." << G4endl;
803     }
804   if ((fSection != "head") && (fSection != "trunk") && (fSection != "full"))
805     G4cout << fSection << " is not defined!" << G4endl;  
806 
807 }
808