ICRP145_HumanPhantoms/src/TETModelImport.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 //
0027 // Author: Haegin Han
0028 // Reference: ICRP Publication 145. Ann. ICRP 49(3), 2020.
0029 // Geant4 Contributors: J. Allison and S. Guatelli
0030 //
0031
0032 #include "TETModelImport.hh"
0033
0034 TETModelImport::TETModelImport(G4bool isAF, G4UIExecutive* ui)
0035 {
0036   // set path for phantom data
0037   char* pPATH = std::getenv("PHANTOM_PATH");
0038   if( pPATH == nullptr )
0039   {
0040     // exception for the case when PHANTOM_PATH environment variable was not set
0041     G4Exception("TETModelImport::TETModelImport","",JustWarning,
0042                 G4String("PHANTOM_PATH environment variable was not set. Using ./ICRP145data.").c_str());
0043     // default path for phantom data
0044     fPhantomDataPath = "./ICRP145data";
0045   }
0046   else
0047   {
0048     // set path for phantom data as PHANTOM_PATH
0049     fPhantomDataPath = pPATH;
0050   }
0051
0052   // set phantom name
0053   if(!isAF) fPhantomName = "MRCP_AM";
0054   else      fPhantomName = "MRCP_AF";
0055
0056   G4cout << "================================================================================"<<G4endl;
0057   G4cout << "\t" << fPhantomName << " was implemented in this CODE!!   "<< G4endl;
0058   G4cout << "================================================================================"<<G4endl;
0059
0060   G4String eleFile      =  fPhantomName + ".ele";
0061   G4String nodeFile     =  fPhantomName + ".node";
0062   G4String materialFile =  fPhantomName + ".material";
0063
0064   // read phantom data files (*. ele, *.node)
0065   DataRead(eleFile, nodeFile);
0066   // read material file (*.material)
0067   MaterialRead(materialFile);
0068   // read colour data file (colour.dat) if this is interactive mode
0069   if(ui) ColourRead();
0070   // print the summary of phantom information
0071   PrintMaterialInfomation();
0072 }
0073
0074 void TETModelImport::DataRead(G4String eleFile, G4String nodeFile)
0075 {
0076   G4String tempStr;
0077   G4int tempInt;
0078
0079   // Read *.node file
0080   //
0081   std::ifstream ifpNode;
0082
0083   ifpNode.open((fPhantomDataPath + "/" + nodeFile).c_str());
0084   if(!ifpNode.is_open())
0085   {
0086     // exception for the case when there is no *.node file
0087     G4Exception("TETModelImport::DataRead","",FatalErrorInArgument,
0088                 G4String("      There is no " + nodeFile ).c_str());
0089   }
0090   G4cout << "  Opening TETGEN node (vertex points: x y z) file '"
0091          << nodeFile << "'" <<G4endl;
0092
0093   G4int numVertex;
0094   G4double xPos, yPos, zPos;
0095   G4double xMin(DBL_MAX), yMin(DBL_MAX), zMin(DBL_MAX);
0096   G4double xMax(DBL_MIN), yMax(DBL_MIN), zMax(DBL_MIN);
0097
0098   ifpNode >> numVertex >> tempInt >> tempInt >> tempInt;
0099
0100   for(G4int i=0; i<numVertex; ++i)
0101   {
0102     ifpNode >> tempInt >> xPos >> yPos >> zPos;
0103
0104     // set the unit
0105     xPos*=cm;
0106     yPos*=cm;
0107     zPos*=cm;
0108
0109     // save the node data as the form of std::vector<G4ThreeVector>
0110     fVertexVector.push_back(G4ThreeVector(xPos, yPos, zPos));
0111
0112     // to get the information of the bounding box of phantom
0113     if (xPos < xMin) xMin = xPos;
0114     if (xPos > xMax) xMax = xPos;
0115     if (yPos < yMin) yMin = yPos;
0116     if (yPos > yMax) yMax = yPos;
0117     if (zPos < zMin) zMin = zPos;
0118     if (zPos > zMax) zMax = zPos;
0119   }
0120
0121   // set the variables for the bounding box and phantom size
0122   fBoundingBox_Min = G4ThreeVector(xMin,yMin,zMin);
0123   fBoundingBox_Max = G4ThreeVector(xMax,yMax,zMax);
0124   fPhantomSize = G4ThreeVector(xMax-xMin,yMax-yMin,zMax-zMin);
0125
0126   ifpNode.close();
0127
0128   // Read *.ele file
0129   //
0130   std::ifstream ifpEle;
0131
0132   ifpEle.open((fPhantomDataPath + "/" + eleFile).c_str());
0133   if(!ifpEle.is_open())
0134   {
0135     // exception for the case when there is no *.ele file
0136     G4Exception("TETModelImport::DataRead","",FatalErrorInArgument,
0137                 G4String("      There is no " + eleFile ).c_str());
0138   }
0139   G4cout << "  Opening TETGEN elements (tetrahedron with node No.) file '"
0140          << eleFile << "'" <<G4endl;
0141
0142   G4int numEle;
0143   ifpEle >> numEle  >> tempInt >> tempInt;
0144
0145   for(G4int i=0; i<numEle; ++i)
0146   {
0147     ifpEle >> tempInt;
0148     auto* ele = new G4int[4];
0149     for(G4int j=0;j<4;j++)
0150     {
0151       ifpEle >> tempInt;
0152       ele[j]=tempInt;
0153     }
0154     fEleVector.push_back(ele);
0155     ifpEle >> tempInt;
0156     fMaterialVector.push_back(tempInt);
0157
0158     // save the element (tetrahedron) data as the form of std::vector<G4Tet*>
0159     fTetVector.push_back(new G4Tet("Tet_Solid",
0160                                    fVertexVector[ele[0]],
0161                                    fVertexVector[ele[1]],
0162                                    fVertexVector[ele[2]],
0163                                    fVertexVector[ele[3]]));
0164
0165     // calculate the total volume and the number of tetrahedrons for each organ
0166     //std::map<G4int, G4double>::iterator FindIter = fVolumeMap.find(fMaterialVector[i]);
0167     auto FindIter = fVolumeMap.find(fMaterialVector[i]);
0168
0169     if(FindIter!=fVolumeMap.end())
0170     {
0171       FindIter->second += fTetVector[i]->GetCubicVolume();
0172       fNumTetMap[fMaterialVector[i]]++;
0173     }
0174     else
0175     {
0176       fVolumeMap[fMaterialVector[i]] = fTetVector[i]->GetCubicVolume();
0177       fNumTetMap[fMaterialVector[i]] = 1;
0178     }
0179   }
0180   ifpEle.close();
0181 }
0182
0183 void TETModelImport::MaterialRead(G4String materialFile)
0184 {
0185   // Read material file (*.material)
0186   //
0187   std::ifstream ifpMat;
0188
0189   ifpMat.open((fPhantomDataPath + "/" + materialFile).c_str());
0190   if(!ifpMat.is_open())
0191   {
0192     // exception for the case when there is no *.material file
0193     G4Exception("TETModelImport::DataRead","",FatalErrorInArgument,
0194     G4String("      There is no " + materialFile ).c_str());
0195   }
0196
0197   G4cout << "  Opening material file '" << materialFile << "'" <<G4endl;
0198
0199   char read_data[50];
0200   char* token;
0201   G4double zaid;
0202   G4double fraction;
0203   G4String MaterialName;
0204   G4double density;
0205   G4int i=0;
0206
0207   while(!ifpMat.eof())
0208   {
0209     ifpMat >> read_data;                   //ex) 'C' RBM
0210     ifpMat >> MaterialName;                //ex)  C 'RBM'
0211     ifpMat >> read_data;
0212     density = std::atof(read_data);        //ex) 1.30
0213     ifpMat >> read_data;                   //ex) g/cm3
0214     ifpMat >> read_data;
0215     token = std::strtok(read_data,"m");
0216     G4int matID = std::atoi(token);        //ex) m'10'
0217     fMaterialIndex.push_back(matID);
0218     fOrganNameMap[matID]= MaterialName;
0219     fDensityMap[matID] = density*g/cm3;
0220
0221     for(i=0 ;  ; ++i)
0222     {
0223       ifpMat >> read_data;
0224       if(std::strcmp(read_data, "C")==0 || ifpMat.eof()) break;
0225
0226       zaid = (G4int)(std::atoi(read_data)/1000);
0227       ifpMat >> read_data;
0228       fraction = -1.0 * std::atof(read_data);
0229       fMaterialIndexMap[matID].push_back(std::make_pair(G4int(zaid), fraction));
0230     }
0231   }
0232   ifpMat.close();
0233
0234   // Construct materials for each organ
0235   //
0236   G4NistManager* nistManager = G4NistManager::Instance();
0237
0238   for(i=0;i<(G4int)fMaterialIndex.size();++i)
0239   {
0240     G4int idx = fMaterialIndex[i];
0241     auto* mat = new G4Material(fOrganNameMap[idx], fDensityMap[idx],
0242                                G4int(fMaterialIndexMap[idx].size()),
0243                                kStateSolid, NTP_Temperature, STP_Pressure);
0244     for(G4int j=0;j<G4int(fMaterialIndexMap[idx].size());++j)
0245     {
0246       mat->AddElement(nistManager->FindOrBuildElement(fMaterialIndexMap[idx][j].first),
0247                                                       fMaterialIndexMap[idx][j].second);
0248     }
0249     fMaterialMap[idx]=mat;
0250     fMassMap[idx]=fDensityMap[idx]*fVolumeMap[idx];
0251   }
0252 }
0253
0254 void TETModelImport::ColourRead()
0255 {
0256   // Read colour data file (colour.dat)
0257   //
0258   std::ifstream ifpColour;
0259
0260   ifpColour.open((fPhantomDataPath + "/" + "colour.dat").c_str());
0261   if(!ifpColour.is_open())
0262   {
0263     // exception for the case when there is no colour.dat file
0264     G4Exception("TETModelImport::DataRead","",FatalErrorInArgument,
0265                 G4String("Colour data file was not found ").c_str());
0266   }
0267
0268   G4cout << "  Opening colour data file 'colour.dat'" <<G4endl;
0269
0270   G4int organID;
0271   G4double red, green, blue, alpha;
0272   while( ifpColour >> organID >> red >> green >> blue >> alpha )
0273   {
0274     fColourMap[organID] = G4Colour(red, green, blue, alpha);
0275   }
0276
0277   ifpColour.close();
0278 }
0279
0280 void TETModelImport::PrintMaterialInfomation()
0281 {
0282   // Print the overall information for each organ
0283   //
0284   G4cout << G4endl
0285          << std::setw(9)  << "Organ ID"
0286          << std::setw(11) << "# of Tet"
0287          << std::setw(11) << "vol [cm3]"
0288          << std::setw(11) << "d [g/cm3]"
0289          << std::setw(11) << "mass [g]"
0290          << "\t" << "organ/tissue"<< G4endl ;
0291
0292   G4cout << "--------------------------------------------------------------------------------"<<G4endl;
0293
0294   std::map<G4int, G4Material*>::const_iterator matIter;
0295   G4cout<<std::setiosflags(std::ios::fixed);
0296   G4cout.precision(3);
0297   for(matIter=fMaterialMap.cbegin(); matIter!=fMaterialMap.cend(); ++matIter)
0298   {
0299     G4int idx = matIter->first;
0300     G4cout << std::setw(9)  << idx                          // organ ID
0301            << std::setw(11) << fNumTetMap[idx]              // # of tetrahedrons
0302            << std::setw(11) << fVolumeMap[idx]/cm3          // organ volume
0303            << std::setw(11) << fMaterialMap[idx] ->GetDensity()/(g/cm3)     // organ density
0304            << std::setw(11) << fMassMap[idx]/g              // organ mass
0305            << "\t"<<fMaterialMap[idx]->GetName() << G4endl; // organ name
0306   }
0307 }