ICRP110_HumanPhantoms/src/ICRP110PhantomConstruction.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // Code developed by:
0027 // S.Guatelli, M. Large and A. Malaroda, University of Wollongong
0028 //
0029 #include "ICRP110PhantomConstruction.hh"
0030 #include "ICRP110PhantomNestedParameterisation.hh"
0031 #include "ICRP110PhantomMaterial_Female.hh"
0032 #include "ICRP110PhantomMaterial_Male.hh"
0033 #include "globals.hh"
0034 #include "G4SystemOfUnits.hh"
0035 #include "G4Box.hh"
0036 #include "G4Colour.hh"
0037 #include "G4LogicalVolume.hh"
0038 #include "G4PVPlacement.hh"
0039 #include "G4PVParameterised.hh"
0040 #include "G4VPhysicalVolume.hh"
0041 #include "G4PVPlacement.hh"
0042 #include "G4PVParameterised.hh"
0043 #include "G4RunManager.hh"
0044 #include "G4VisAttributes.hh"
0045 #include <map>
0046 #include <cstdlib>
0047
0048 ICRP110PhantomConstruction::ICRP110PhantomConstruction():
0049    fMotherVolume(nullptr), fPhantomContainer(nullptr),
0050    fNVoxelX(0), fNVoxelY(0), fNVoxelZ(0),
0051    fVoxelHalfDimX(0), fVoxelHalfDimY(0), fVoxelHalfDimZ(0),
0052    fMinX(0),fMaxX(0), fMinY(0), fMaxY(0),
0053    fMinZ(0), fMaxZ(0), fNoFiles(0), fNVoxels(0),
0054    fMateIDs(nullptr)
0055 {
0056   fMessenger = new ICRP110PhantomMessenger(this);
0057   // the messenger allows to set the sex of the phantom
0058   // interactively
0059   fMaterial_Female = new ICRP110PhantomMaterial_Female();
0060   fMaterial_Male = new ICRP110PhantomMaterial_Male();
0061   fSex = "female"; // Female phantom is the default option
0062   fSection = "head"; // Head partial phantom is the default option
0063 }
0064
0065 ICRP110PhantomConstruction::~ICRP110PhantomConstruction()
0066 {
0067   delete fMaterial_Female;
0068   delete fMaterial_Male;
0069   delete fMessenger;
0070 }
0071
0072 G4VPhysicalVolume* ICRP110PhantomConstruction::Construct()
0073 {
0074   // Define Material Air
0075   G4double A;  // atomic mass
0076   G4double Z;  // atomic number
0077   G4double d;  // density
0078
0079   A = 14.01*g/mole;
0080   auto elN = new G4Element("Nitrogen","N",Z = 7.,A);
0081   A = 16.00*g/mole;
0082   auto elO = new G4Element("Oxygen","O",Z = 8.,A);
0083
0084   d = 0.001 *g/cm3;
0085   auto matAir = new G4Material("Air",d,2);
0086   matAir -> AddElement(elN,0.8);
0087   matAir -> AddElement(elO,0.2);
0088
0089  std::vector<G4Material*> pMaterials;
0090
0091  if(fSex == "female"){
0092
0093   fMaterial_Female -> DefineMaterials();
0094 //----- Store materials in a vector
0095     pMaterials.push_back(matAir);
0096     pMaterials.push_back(fMaterial_Female -> GetMaterial("teeth"));
0097     pMaterials.push_back(fMaterial_Female -> GetMaterial("bone"));
0098     pMaterials.push_back(fMaterial_Female -> GetMaterial("humeri_upper"));
0099     pMaterials.push_back(fMaterial_Female -> GetMaterial("humeri_lower"));
0100     pMaterials.push_back(fMaterial_Female -> GetMaterial("arm_lower"));
0101     pMaterials.push_back(fMaterial_Female -> GetMaterial("hand"));
0102     pMaterials.push_back(fMaterial_Female -> GetMaterial("clavicle"));
0103     pMaterials.push_back(fMaterial_Female -> GetMaterial("cranium"));
0104     pMaterials.push_back(fMaterial_Female -> GetMaterial("femora_upper"));
0105     pMaterials.push_back(fMaterial_Female -> GetMaterial("femora_lower"));
0106     pMaterials.push_back(fMaterial_Female -> GetMaterial("leg_lower"));
0107     pMaterials.push_back(fMaterial_Female -> GetMaterial("foot"));
0108     pMaterials.push_back(fMaterial_Female -> GetMaterial("mandible"));
0109     pMaterials.push_back(fMaterial_Female -> GetMaterial("pelvis"));
0110     pMaterials.push_back(fMaterial_Female -> GetMaterial("ribs"));
0111     pMaterials.push_back(fMaterial_Female -> GetMaterial("scapulae"));
0112     pMaterials.push_back(fMaterial_Female -> GetMaterial("spine_cervical"));
0113     pMaterials.push_back(fMaterial_Female -> GetMaterial("spine_thoratic"));
0114     pMaterials.push_back(fMaterial_Female -> GetMaterial("spine_lumbar"));
0115     pMaterials.push_back(fMaterial_Female -> GetMaterial("sacrum"));
0116     pMaterials.push_back(fMaterial_Female -> GetMaterial("sternum"));
0117     pMaterials.push_back(fMaterial_Female -> GetMaterial("hf_upper"));
0118     pMaterials.push_back(fMaterial_Female -> GetMaterial("hf_lower"));
0119     pMaterials.push_back(fMaterial_Female -> GetMaterial("med_lowerarm"));
0120     pMaterials.push_back(fMaterial_Female -> GetMaterial("med_lowerleg"));
0121     pMaterials.push_back(fMaterial_Female -> GetMaterial("cartilage"));
0122     pMaterials.push_back(fMaterial_Female -> GetMaterial("skin"));
0123     pMaterials.push_back(fMaterial_Female -> GetMaterial("blood"));
0124     pMaterials.push_back(fMaterial_Female -> GetMaterial("muscle"));
0125     pMaterials.push_back(fMaterial_Female -> GetMaterial("liver"));
0126     pMaterials.push_back(fMaterial_Female -> GetMaterial("pancreas"));
0127     pMaterials.push_back(fMaterial_Female -> GetMaterial("brain"));
0128     pMaterials.push_back(fMaterial_Female -> GetMaterial("heart"));
0129     pMaterials.push_back(fMaterial_Female -> GetMaterial("eye"));
0130     pMaterials.push_back(fMaterial_Female -> GetMaterial("kidney"));
0131     pMaterials.push_back(fMaterial_Female -> GetMaterial("stomach"));
0132     pMaterials.push_back(fMaterial_Female -> GetMaterial("intestine_sml"));
0133     pMaterials.push_back(fMaterial_Female -> GetMaterial("intestine_lrg"));
0134     pMaterials.push_back(fMaterial_Female -> GetMaterial("spleen"));
0135     pMaterials.push_back(fMaterial_Female -> GetMaterial("thyroid"));
0136     pMaterials.push_back(fMaterial_Female -> GetMaterial("bladder"));
0137     pMaterials.push_back(fMaterial_Female -> GetMaterial("ovaries_testes"));
0138     pMaterials.push_back(fMaterial_Female -> GetMaterial("adrenals"));
0139     pMaterials.push_back(fMaterial_Female -> GetMaterial("oesophagus"));
0140     pMaterials.push_back(fMaterial_Female -> GetMaterial("misc"));
0141     pMaterials.push_back(fMaterial_Female -> GetMaterial("uterus_prostate"));
0142     pMaterials.push_back(fMaterial_Female -> GetMaterial("lymph"));
0143     pMaterials.push_back(fMaterial_Female -> GetMaterial("breast_glandular"));
0144     pMaterials.push_back(fMaterial_Female -> GetMaterial("breast_adipose"));
0145     pMaterials.push_back(fMaterial_Female -> GetMaterial("lung"));
0146     pMaterials.push_back(fMaterial_Female -> GetMaterial("gastro_content"));
0147     pMaterials.push_back(fMaterial_Female -> GetMaterial("urine"));
0148  }
0149  else if (fSex == "male"){
0150  // MATT do the same here
0151     fMaterial_Male -> DefineMaterials();
0152
0153 //----- Store materials in a vector
0154     pMaterials.push_back(matAir);
0155     pMaterials.push_back(fMaterial_Male -> GetMaterial("teeth"));
0156     pMaterials.push_back(fMaterial_Male -> GetMaterial("bone"));
0157     pMaterials.push_back(fMaterial_Male -> GetMaterial("humeri_upper"));
0158     pMaterials.push_back(fMaterial_Male -> GetMaterial("humeri_lower"));
0159     pMaterials.push_back(fMaterial_Male -> GetMaterial("arm_lower"));
0160     pMaterials.push_back(fMaterial_Male -> GetMaterial("hand"));
0161     pMaterials.push_back(fMaterial_Male -> GetMaterial("clavicle"));
0162     pMaterials.push_back(fMaterial_Male -> GetMaterial("cranium"));
0163     pMaterials.push_back(fMaterial_Male -> GetMaterial("femora_upper"));
0164     pMaterials.push_back(fMaterial_Male -> GetMaterial("femora_lower"));
0165     pMaterials.push_back(fMaterial_Male -> GetMaterial("leg_lower"));
0166     pMaterials.push_back(fMaterial_Male -> GetMaterial("foot"));
0167     pMaterials.push_back(fMaterial_Male -> GetMaterial("mandible"));
0168     pMaterials.push_back(fMaterial_Male -> GetMaterial("pelvis"));
0169     pMaterials.push_back(fMaterial_Male -> GetMaterial("ribs"));
0170     pMaterials.push_back(fMaterial_Male -> GetMaterial("scapulae"));
0171     pMaterials.push_back(fMaterial_Male -> GetMaterial("spine_cervical"));
0172     pMaterials.push_back(fMaterial_Male -> GetMaterial("spine_thoratic"));
0173     pMaterials.push_back(fMaterial_Male -> GetMaterial("spine_lumbar"));
0174     pMaterials.push_back(fMaterial_Male -> GetMaterial("sacrum"));
0175     pMaterials.push_back(fMaterial_Male -> GetMaterial("sternum"));
0176     pMaterials.push_back(fMaterial_Male -> GetMaterial("hf_upper"));
0177     pMaterials.push_back(fMaterial_Male -> GetMaterial("hf_lower"));
0178     pMaterials.push_back(fMaterial_Male -> GetMaterial("med_lowerarm"));
0179     pMaterials.push_back(fMaterial_Male -> GetMaterial("med_lowerleg"));
0180     pMaterials.push_back(fMaterial_Male -> GetMaterial("cartilage"));
0181     pMaterials.push_back(fMaterial_Male -> GetMaterial("skin"));
0182     pMaterials.push_back(fMaterial_Male -> GetMaterial("blood"));
0183     pMaterials.push_back(fMaterial_Male -> GetMaterial("muscle"));
0184     pMaterials.push_back(fMaterial_Male -> GetMaterial("liver"));
0185     pMaterials.push_back(fMaterial_Male -> GetMaterial("pancreas"));
0186     pMaterials.push_back(fMaterial_Male -> GetMaterial("brain"));
0187     pMaterials.push_back(fMaterial_Male -> GetMaterial("heart"));
0188     pMaterials.push_back(fMaterial_Male -> GetMaterial("eye"));
0189     pMaterials.push_back(fMaterial_Male -> GetMaterial("kidney"));
0190     pMaterials.push_back(fMaterial_Male -> GetMaterial("stomach"));
0191     pMaterials.push_back(fMaterial_Male -> GetMaterial("intestine_sml"));
0192     pMaterials.push_back(fMaterial_Male -> GetMaterial("intestine_lrg"));
0193     pMaterials.push_back(fMaterial_Male -> GetMaterial("spleen"));
0194     pMaterials.push_back(fMaterial_Male -> GetMaterial("thyroid"));
0195     pMaterials.push_back(fMaterial_Male -> GetMaterial("bladder"));
0196     pMaterials.push_back(fMaterial_Male -> GetMaterial("ovaries_testes"));
0197     pMaterials.push_back(fMaterial_Male -> GetMaterial("adrenals"));
0198     pMaterials.push_back(fMaterial_Male -> GetMaterial("oesophagus"));
0199     pMaterials.push_back(fMaterial_Male -> GetMaterial("misc"));
0200     pMaterials.push_back(fMaterial_Male -> GetMaterial("uterus_prostate"));
0201     pMaterials.push_back(fMaterial_Male -> GetMaterial("lymph"));
0202     pMaterials.push_back(fMaterial_Male -> GetMaterial("breast_glandular"));
0203     pMaterials.push_back(fMaterial_Male -> GetMaterial("breast_adipose"));
0204     pMaterials.push_back(fMaterial_Male -> GetMaterial("lung"));
0205     pMaterials.push_back(fMaterial_Male -> GetMaterial("gastro_content"));
0206     pMaterials.push_back(fMaterial_Male -> GetMaterial("urine"));
0207
0208 }
0209
0210   // World Volume
0211   G4double worldSize = 2.*m ;
0212   G4Box* world = new G4Box("world", worldSize, worldSize, worldSize);
0213
0214   auto logicWorld = new G4LogicalVolume(world,
0215                              matAir,
0216                          "logicalWorld", nullptr, nullptr,nullptr);
0217
0218   fMotherVolume = new G4PVPlacement(nullptr,G4ThreeVector(),
0219                     "physicalWorld",
0220                     logicWorld,
0221                     nullptr,
0222                     false,
0223                     0);
0224
0225   logicWorld -> SetVisAttributes(G4VisAttributes::GetInvisible());
0226
0227   G4cout << "World has been built" << G4endl;
0228
0229   G4cout << "Phantom Sex: " << fSex << G4endl;
0230   G4cout << "Phantom Section: " << fSection << G4endl;
0231   ReadPhantomData(fSex, fSection);
0232
0233   G4cout << "Number of X,Y,Z voxels = " << fNVoxelX << ", " << fNVoxelY << ", " << fNVoxelZ << G4endl;
0234
0235 //----- Define the volume that contains all the voxels
0236   G4Box* fContainer_solid = new G4Box("phantomContainer",fNVoxelX*fVoxelHalfDimX*mm,
0237                                fNVoxelY*fVoxelHalfDimY*mm,
0238                                fNVoxelZ*fVoxelHalfDimZ*mm);
0239
0240   auto fContainer_logic = new G4LogicalVolume( fContainer_solid,
0241                                                             matAir,
0242                                                             "phantomContainer",
0243                                                              nullptr, nullptr, nullptr);
0244
0245   fMaxX = fNVoxelX*fVoxelHalfDimX*mm; // Max X along X axis of the voxelised geometry
0246   fMaxY = fNVoxelY*fVoxelHalfDimY*mm; // Max Y
0247   fMaxZ = fNVoxelZ*fVoxelHalfDimZ*mm; // Max Z
0248
0249   fMinX = -fNVoxelX*fVoxelHalfDimX*mm;// Min X
0250   fMinY = -fNVoxelY*fVoxelHalfDimY*mm;// Min Y
0251   fMinZ = -fNVoxelZ*fVoxelHalfDimZ*mm;// Min Z
0252
0253   G4ThreeVector posCentreVoxels((fMinX+fMaxX)/2.,(fMinY+fMaxY)/2.,(fMinZ+fMaxZ)/2.);
0254
0255   G4cout << " placing voxel container volume at " << posCentreVoxels << G4endl;
0256
0257
0258   fPhantomContainer
0259   = new G4PVPlacement(nullptr,                     // rotation
0260                       posCentreVoxels,
0261                       fContainer_logic,     // The logic volume
0262                       "phantomContainer",  // Name
0263                       logicWorld,         // Mother
0264                       false,            // No op. bool.
0265                       1);              // Copy number
0266
0267   fContainer_logic -> SetVisAttributes(new G4VisAttributes(G4Colour(1.,0.,0.,0.)));
0268
0269
0270 // Define the voxelised phantom here
0271 // Replication of air Phantom Volume.
0272
0273 //--- Slice the phantom along Y axis
0274    G4String yRepName("RepY");
0275    G4VSolid* solYRep = new G4Box(yRepName,fNVoxelX*fVoxelHalfDimX,
0276                                   fVoxelHalfDimY, fNVoxelZ*fVoxelHalfDimZ);
0277    auto logYRep = new G4LogicalVolume(solYRep,matAir,yRepName);
0278    new G4PVReplica(yRepName,logYRep,fContainer_logic,kYAxis, fNVoxelY,fVoxelHalfDimY*2.);
0279
0280    logYRep -> SetVisAttributes(new G4VisAttributes(G4VisAttributes::GetInvisible()));
0281
0282 //--- Slice the phantom along X axis
0283    G4String xRepName("RepX");
0284    G4VSolid* solXRep = new G4Box(xRepName,fVoxelHalfDimX,fVoxelHalfDimY,
0285                                   fNVoxelZ*fVoxelHalfDimZ);
0286    auto logXRep = new G4LogicalVolume(solXRep,matAir,xRepName);
0287    new G4PVReplica(xRepName,logXRep,logYRep,kXAxis,fNVoxelX,fVoxelHalfDimX*2.);
0288
0289    logXRep -> SetVisAttributes(new G4VisAttributes(G4VisAttributes::GetInvisible()));
0290
0291    //----- Voxel solid and logical volumes
0292    //--- Slice along Z axis
0293    G4VSolid* solidVoxel = new G4Box("phantom",fVoxelHalfDimX, fVoxelHalfDimY,fVoxelHalfDimZ);
0294    auto logicVoxel = new G4LogicalVolume(solidVoxel,matAir,"phantom");
0295
0296    logicVoxel -> SetVisAttributes(new G4VisAttributes(G4VisAttributes::GetInvisible()));
0297
0298    // Parameterisation to define the material of each voxel
0299    G4ThreeVector halfVoxelSize(fVoxelHalfDimX,fVoxelHalfDimY,fVoxelHalfDimZ);
0300
0301    auto param =  new ICRP110PhantomNestedParameterisation(halfVoxelSize, pMaterials);
0302
0303    new G4PVParameterised("phantom",    // their name
0304                           logicVoxel, // their logical volume
0305                           logXRep,      // Mother logical volume
0306                           kZAxis,       // Are placed along this axis
0307                           fNVoxelZ,      // Number of cells
0308                           param);       // Parameterisation
0309
0310     param -> SetMaterialIndices(fMateIDs); // fMateIDs is  the vector with Material ID associated to each voxel, from ASCII input data files.
0311     param -> SetNoVoxel(fNVoxelX,fNVoxelY,fNVoxelZ);
0312
0313   return fMotherVolume;
0314 }
0315
0316 void ICRP110PhantomConstruction::ReadPhantomData(const G4String& sex, const G4String& section)
0317 {
0318
0319   // This method reads the information of ICRPdata/FemaleData.dat or
0320   // ICRPdata/MaleData.data depending on the sex of the chosen phantom
0321
0322 fSex = sex;
0323 fSection = section;
0324
0325 G4String dataFile;
0326
0327     if (fSex == "female")
0328     {
0329         if (fSection == "head")
0330         {
0331           dataFile = "ICRPdata/FemaleHead.dat";
0332         }
0333         else if (fSection == "trunk")
0334         {
0335           dataFile = "ICRPdata/FemaleTrunk.dat";
0336         }
0337         else if (fSection == "full")
0338         {
0339           dataFile = "ICRPdata/FemaleData.dat";
0340         }
0341     }
0342     if (fSex == "male")
0343     {
0344         if (fSection == "head")
0345         {
0346           dataFile = "ICRPdata/MaleHead.dat";
0347         }
0348         else if (fSection == "trunk")
0349         {
0350           dataFile = "ICRPdata/MaleTrunk.dat";
0351         }
0352         else if (fSection == "full")
0353         {
0354           dataFile = "ICRPdata/MaleData.dat";
0355         }
0356     }
0357
0358     G4cout << "Data file " << dataFile << " is read by Detector Construction." << G4endl;
0359
0360 // The data.dat file in directory/build/ICRPdata/ contains the information
0361 // to build the phantoms. For more details look in the README file.
0362
0363   //input file named finDF which consists of dataFile as a string object
0364   std::ifstream finDF(dataFile.c_str());
0365
0366
0367   G4String fname;
0368
0369 if(finDF.good() != 1 ) //check that the file is good and working
0370  {
0371   G4String descript = "Problem reading data file: "+dataFile;
0372   G4Exception(" HumanPhantomConstruction::ReadPhantomData"," ",
0373               FatalException,descript);
0374   }
0375
0376   finDF >> fNoFiles;
0377   G4cout << "Number of files = " << fNoFiles << G4endl;
0378   finDF >> fNVoxelX;      //Inputs number of X-Voxels
0379   finDF >> fNVoxelY;      //Y-Voxels
0380   fNVoxelZ = fNoFiles;    //Z-Voxels (equal to number of slice files built/read)
0381   finDF >> fVoxelHalfDimX;
0382   finDF >> fVoxelHalfDimY;
0383   finDF >> fVoxelHalfDimZ;
0384   G4cout << "Number of X,Y,Z voxels = " << fNVoxelX << ", " << fNVoxelY << ", " << fNVoxelZ <<G4endl;
0385
0386   fNVoxels = fNVoxelX*fNVoxelY*fNVoxelZ;
0387   G4cout << "Total Number of Voxels = " << fNVoxels << G4endl;
0388
0389   G4int nMaterials;
0390   finDF >> nMaterials;
0391   G4String mateName;
0392   G4int nmate;
0393
0394 //-----Read materials and associate with material ID number------//
0395
0396   for( G4int ii = 0; ii < nMaterials; ii++ ){
0397     finDF >> nmate;
0398     finDF >> mateName;
0399
0400     // This allows to skip empty spaces and tabs in the string
0401     if( mateName[0] == '"' && mateName[mateName.length()-1] == '"' )
0402     {
0403       mateName = mateName.substr(1,mateName.length()-2);
0404     }
0405
0406     // To uncomment for eventual debugging
0407     /* G4cout << "GmReadPhantomG4Geometry::ReadPhantomData reading nmate "
0408            << ii << " = " << nmate
0409            << " mate " << mateName << G4endl;*/
0410
0411     if( ii != nmate ) {
0412     G4Exception("GmReadPhantomG4Geometry::ReadPhantomData",
0413                 "Wrong argument",
0414                 FatalErrorInArgument,
0415                 "Material number should be in increasing order:wrong material number");
0416                 }
0417       }
0418
0419 fMateIDs = new size_t[fNVoxels]; //Array with Material ID for each voxel
0420
0421 G4cout << "ICRP110PhantomConstruction::ReadPhantomDataFile is openining the following phantom files: " << G4endl;
0422
0423 for(G4int i = 0; i < fNoFiles; i++ )
0424   {
0425     finDF >> fname;
0426     ReadPhantomDataFile(fSex, fname, i);
0427   }
0428
0429 finDF.close();
0430 }
0431
0432 //----------------Opens phantom ASCII slice files to construct the phantom from-----------------//
0433
0434 void ICRP110PhantomConstruction::ReadPhantomDataFile(const G4String& sex, const G4String& fileName, G4int numberFile)
0435 {
0436 G4cout << fileName << G4endl;
0437
0438 fSex = sex;
0439
0440 G4String slice;
0441
0442     if (fSex == "female")
0443     {
0444       slice = "ICRPdata/ICRP110_g4dat/AF/"+fileName;
0445     }
0446     if (fSex == "male")
0447     {
0448       slice = "ICRPdata/ICRP110_g4dat/AM/"+fileName;
0449     }
0450
0451   std::ifstream fin(slice.c_str(), std::ios_base::in);
0452
0453   if( !fin.is_open() ) {
0454     G4Exception("HumanPhantomConstruction::ReadPhantomDataFile",
0455                 "",
0456                 FatalErrorInArgument,
0457                 G4String("File not found " + fileName ).c_str());
0458   }
0459
0460     for( G4int iy = 0; iy < fNVoxelY; iy++ ) {
0461       for( G4int ix = 0; ix < fNVoxelX; ix++ ) {
0462       if (ix == 0 && iy == 0)
0463         {
0464           G4int dudX,dudY,dudZ;
0465           fin >> dudX >> dudY >> dudZ ;
0466           // Dummy method to skip the first three lines of the files
0467           // which are not used here
0468         }
0469
0470         G4int nnew = ix + (iy)*fNVoxelX + numberFile*fNVoxelX*fNVoxelY;
0471         G4int OrgID;
0472         fin >> OrgID;
0473
0474         G4int mateID_out;
0475
0476 // The code below associates organ ID numbers (called here mateID) from ASCII slice
0477 // files with material ID numbers (called here mateID_out) as defined in ICRP110PhantomMaterials
0478 // Material and Organ IDs are associated as stated in AM_organs.dat and FM_organs.dat depending on
0479 // the sex of the phantom (male and female, respctively)
0480
0481     if (OrgID==128)
0482     {
0483     mateID_out=1;
0484     }
0485
0486     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)
0487     {
0488     mateID_out=2;
0489     }
0490
0491     else if (OrgID==14)
0492     {
0493     mateID_out=3;
0494     }
0495
0496     else if (OrgID==17)
0497     {
0498     mateID_out=4;
0499     }
0500
0501     else if (OrgID==20)
0502     {
0503     mateID_out=5;
0504     }
0505
0506     else if (OrgID==23)
0507     {
0508     mateID_out=6;
0509     }
0510
0511     else if (OrgID==25)
0512     {
0513     mateID_out=7;
0514     }
0515
0516     else if (OrgID==27)
0517     {
0518     mateID_out=8;
0519     }
0520
0521     else if (OrgID==29)
0522     {
0523     mateID_out=9;
0524     }
0525
0526     else if (OrgID==32)
0527     {
0528     mateID_out=10;
0529     }
0530
0531     else if (OrgID==35)
0532     {
0533     mateID_out=11;
0534     }
0535
0536     else if (OrgID==38)
0537     {
0538     mateID_out=12;
0539     }
0540
0541     else if (OrgID==40)
0542     {
0543     mateID_out=13;
0544     }
0545
0546     else if (OrgID==42)
0547     {
0548     mateID_out=14;
0549     }
0550
0551     else if (OrgID==44)
0552     {
0553     mateID_out=15;
0554     }
0555
0556     else if (OrgID==46)
0557     {
0558     mateID_out=16;
0559     }
0560
0561     else if (OrgID==48)
0562     {
0563     mateID_out=17;
0564     }
0565
0566     else if (OrgID==50)
0567     {
0568     mateID_out=18;
0569     }
0570
0571     else if (OrgID==52)
0572     {
0573     mateID_out=19;
0574     }
0575
0576     else if (OrgID==54)
0577     {
0578     mateID_out=20;
0579     }
0580
0581     else if (OrgID==56)
0582     {
0583     mateID_out=21;
0584     }
0585
0586     else if (OrgID==15 || OrgID==30)
0587     {
0588     mateID_out=22;
0589     }
0590
0591     else if (OrgID==18 || OrgID==33)
0592     {
0593     mateID_out=23;
0594     }
0595
0596     else if (OrgID==21)
0597     {
0598     mateID_out=24;
0599     }
0600
0601     else if (OrgID==36)
0602     {
0603     mateID_out=25;
0604     }
0605
0606     else if (OrgID==57 || OrgID==58 || OrgID==59 || OrgID==60)
0607     {
0608     mateID_out=26;
0609     }
0610
0611     else if (OrgID==122 || OrgID==123 || OrgID==124 || OrgID==125 || OrgID==141 )
0612     {
0613     mateID_out=27;
0614     }
0615
0616     else if (OrgID==9 || OrgID==10 || OrgID==11 || OrgID==12 || OrgID==88 || OrgID==96 || OrgID==98)
0617     {
0618     mateID_out=28;
0619     }
0620
0621     else if (OrgID==5 || OrgID==6 || OrgID==106 || OrgID==107 || OrgID==108 || OrgID==109 || OrgID==133)
0622     {
0623     mateID_out=29;
0624     }
0625
0626     else if (OrgID==95)
0627     {
0628     mateID_out=30;
0629     }
0630
0631     else if (OrgID==113)
0632     {
0633     mateID_out=31;
0634     }
0635
0636     else if (OrgID==61)
0637     {
0638     mateID_out=32;
0639     }
0640
0641     else if (OrgID==87)
0642     {
0643     mateID_out=33;
0644     }
0645
0646     else if (OrgID==66 || OrgID==67 || OrgID==68 || OrgID==69)
0647     {
0648     mateID_out=34;
0649     }
0650
0651     else if (OrgID==89 || OrgID==90 || OrgID==91 || OrgID==92 || OrgID==93 || OrgID==94)
0652     {
0653     mateID_out=35;
0654     }
0655
0656     else if (OrgID==72)
0657     {
0658     mateID_out=36;
0659     }
0660
0661     else if (OrgID==74)
0662     {
0663     mateID_out=37;
0664     }
0665
0666     else if (OrgID==76 || OrgID==78 || OrgID==80 || OrgID==82 || OrgID==84 || OrgID==86)
0667     {
0668     mateID_out=38;
0669     }
0670
0671     else if (OrgID==127)
0672     {
0673     mateID_out=39;
0674     }
0675
0676     else if (OrgID==132)
0677     {
0678     mateID_out=40;
0679     }
0680
0681     else if (OrgID==137)
0682     {
0683     mateID_out=41;
0684     }
0685
0686     else if (OrgID==111 || OrgID==112 || OrgID==129 || OrgID==130)
0687     {
0688     mateID_out=42;
0689     }
0690
0691     else if (OrgID==1 || OrgID==2)
0692     {
0693     mateID_out=43;
0694     }
0695
0696     else if (OrgID==110)
0697     {
0698     mateID_out=44;
0699     }
0700
0701     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)
0702   {
0703   mateID_out=45;
0704   }
0705
0706     else if (OrgID==115 || OrgID==139)
0707     {
0708     mateID_out=46;
0709     }
0710
0711     else if (OrgID==100 || OrgID==101 || OrgID==102 || OrgID==103 || OrgID==104 || OrgID==105)
0712     {
0713     mateID_out=47;
0714     }
0715
0716     else if (OrgID==63 || OrgID==65)
0717     {
0718     mateID_out=48;
0719     }
0720
0721     else if (OrgID==62 || OrgID==64 || OrgID==116 || OrgID==117 || OrgID==118 || OrgID==119)
0722     {
0723     mateID_out=49;
0724     }
0725
0726     else if (OrgID==97 || OrgID==99)
0727     {
0728     mateID_out=50;
0729     }
0730
0731     else if (OrgID==73 || OrgID==75 || OrgID==77 || OrgID==79 || OrgID==81 || OrgID==83 || OrgID==85)
0732     {
0733     mateID_out=51;
0734     }
0735
0736     else if (OrgID==138)
0737     {
0738     mateID_out=52;
0739     }
0740
0741     else if (OrgID==0 || OrgID==140)
0742     {
0743     mateID_out=0;
0744     }
0745
0746     else
0747     {
0748     mateID_out=OrgID;
0749     }
0750
0751         G4int nMaterials = 53;
0752         if( mateID_out < 0 || mateID_out >= nMaterials ) {
0753           G4Exception("GmReadPhantomG4Geometry::ReadPhantomData",
0754                       "Wrong index in phantom file",
0755                       FatalException,
0756                       G4String("It should be between 0 and "
0757                               + G4UIcommand::ConvertToString(nMaterials-1)
0758                               + ", while it is "
0759                               + G4UIcommand::ConvertToString(OrgID)).c_str());
0760
0761 //-------------Store Material IDs and position/reference number within phantom in vector---------------//
0762      }
0763
0764           fMateIDs[nnew] = mateID_out;
0765
0766
0767       }
0768    }
0769 }
0770
0771 //-----------Define phantom sex (male or female) to be constructed-------------//
0772 void ICRP110PhantomConstruction::SetPhantomSex(G4String newSex)
0773 {
0774   fSex = newSex;
0775
0776   if (fSex == "male")
0777     {
0778       G4cout << ">> Male Phantom will be built." << G4endl;
0779     }
0780   if (fSex == "female")
0781     {
0782       G4cout << ">> Female Phantom will be built." << G4endl;
0783     }
0784   if ((fSex != "female") && (fSex != "male"))
0785     G4cout << fSex << " is not defined!" << G4endl;
0786 }
0787
0788 //-----------Define phantom section to be constructed-------------//
0789 void ICRP110PhantomConstruction::SetPhantomSection(G4String newSection)
0790 {
0791   fSection = newSection;
0792   if (fSection == "head")
0793     {
0794       G4cout << ">> Partial Head Phantom will be built." << G4endl;
0795     }
0796   if (fSection == "trunk")
0797     {
0798       G4cout << ">> Partial Trunk Phantom will be built." << G4endl;
0799     }
0800   if (fSection == "full")
0801     {
0802       G4cout << ">> Full/Custom Phantom will be built." << G4endl;
0803     }
0804   if ((fSection != "head") && (fSection != "trunk") && (fSection != "full"))
0805     G4cout << fSection << " is not defined!" << G4endl;
0806
0807 }