| Geant4 Cross Reference |
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 //
27 // Author: Haegin Han
28 // Reference: ICRP Publication 145. Ann. ICRP 49(3), 2020.
29 // Geant4 Contributors: J. Allison and S. Guatelli
30 //
31
32 #include "TETModelImport.hh"
33
34 TETModelImport::TETModelImport(G4bool isAF, G4UIExecutive* ui)
35 {
36 // set path for phantom data
37 char* pPATH = std::getenv("PHANTOM_PATH");
38 if( pPATH == nullptr )
39 {
40 // exception for the case when PHANTOM_PATH environment variable was not set
41 G4Exception("TETModelImport::TETModelImport","",JustWarning,
42 G4String("PHANTOM_PATH environment variable was not set. Using ./ICRP145data.").c_str());
43 // default path for phantom data
44 fPhantomDataPath = "./ICRP145data";
45 }
46 else
47 {
48 // set path for phantom data as PHANTOM_PATH
49 fPhantomDataPath = pPATH;
50 }
51
52 // set phantom name
53 if(!isAF) fPhantomName = "MRCP_AM";
54 else fPhantomName = "MRCP_AF";
55
56 G4cout << "================================================================================"<<G4endl;
57 G4cout << "\t" << fPhantomName << " was implemented in this CODE!! "<< G4endl;
58 G4cout << "================================================================================"<<G4endl;
59
60 G4String eleFile = fPhantomName + ".ele";
61 G4String nodeFile = fPhantomName + ".node";
62 G4String materialFile = fPhantomName + ".material";
63
64 // read phantom data files (*. ele, *.node)
65 DataRead(eleFile, nodeFile);
66 // read material file (*.material)
67 MaterialRead(materialFile);
68 // read colour data file (colour.dat) if this is interactive mode
69 if(ui) ColourRead();
70 // print the summary of phantom information
71 PrintMaterialInfomation();
72 }
73
74 void TETModelImport::DataRead(G4String eleFile, G4String nodeFile)
75 {
76 G4String tempStr;
77 G4int tempInt;
78
79 // Read *.node file
80 //
81 std::ifstream ifpNode;
82
83 ifpNode.open((fPhantomDataPath + "/" + nodeFile).c_str());
84 if(!ifpNode.is_open())
85 {
86 // exception for the case when there is no *.node file
87 G4Exception("TETModelImport::DataRead","",FatalErrorInArgument,
88 G4String(" There is no " + nodeFile ).c_str());
89 }
90 G4cout << " Opening TETGEN node (vertex points: x y z) file '"
91 << nodeFile << "'" <<G4endl;
92
93 G4int numVertex;
94 G4double xPos, yPos, zPos;
95 G4double xMin(DBL_MAX), yMin(DBL_MAX), zMin(DBL_MAX);
96 G4double xMax(DBL_MIN), yMax(DBL_MIN), zMax(DBL_MIN);
97
98 ifpNode >> numVertex >> tempInt >> tempInt >> tempInt;
99
100 for(G4int i=0; i<numVertex; ++i)
101 {
102 ifpNode >> tempInt >> xPos >> yPos >> zPos;
103
104 // set the unit
105 xPos*=cm;
106 yPos*=cm;
107 zPos*=cm;
108
109 // save the node data as the form of std::vector<G4ThreeVector>
110 fVertexVector.push_back(G4ThreeVector(xPos, yPos, zPos));
111
112 // to get the information of the bounding box of phantom
113 if (xPos < xMin) xMin = xPos;
114 if (xPos > xMax) xMax = xPos;
115 if (yPos < yMin) yMin = yPos;
116 if (yPos > yMax) yMax = yPos;
117 if (zPos < zMin) zMin = zPos;
118 if (zPos > zMax) zMax = zPos;
119 }
120
121 // set the variables for the bounding box and phantom size
122 fBoundingBox_Min = G4ThreeVector(xMin,yMin,zMin);
123 fBoundingBox_Max = G4ThreeVector(xMax,yMax,zMax);
124 fPhantomSize = G4ThreeVector(xMax-xMin,yMax-yMin,zMax-zMin);
125
126 ifpNode.close();
127
128 // Read *.ele file
129 //
130 std::ifstream ifpEle;
131
132 ifpEle.open((fPhantomDataPath + "/" + eleFile).c_str());
133 if(!ifpEle.is_open())
134 {
135 // exception for the case when there is no *.ele file
136 G4Exception("TETModelImport::DataRead","",FatalErrorInArgument,
137 G4String(" There is no " + eleFile ).c_str());
138 }
139 G4cout << " Opening TETGEN elements (tetrahedron with node No.) file '"
140 << eleFile << "'" <<G4endl;
141
142 G4int numEle;
143 ifpEle >> numEle >> tempInt >> tempInt;
144
145 for(G4int i=0; i<numEle; ++i)
146 {
147 ifpEle >> tempInt;
148 auto* ele = new G4int[4];
149 for(G4int j=0;j<4;j++)
150 {
151 ifpEle >> tempInt;
152 ele[j]=tempInt;
153 }
154 fEleVector.push_back(ele);
155 ifpEle >> tempInt;
156 fMaterialVector.push_back(tempInt);
157
158 // save the element (tetrahedron) data as the form of std::vector<G4Tet*>
159 fTetVector.push_back(new G4Tet("Tet_Solid",
160 fVertexVector[ele[0]],
161 fVertexVector[ele[1]],
162 fVertexVector[ele[2]],
163 fVertexVector[ele[3]]));
164
165 // calculate the total volume and the number of tetrahedrons for each organ
166 //std::map<G4int, G4double>::iterator FindIter = fVolumeMap.find(fMaterialVector[i]);
167 auto FindIter = fVolumeMap.find(fMaterialVector[i]);
168
169 if(FindIter!=fVolumeMap.end())
170 {
171 FindIter->second += fTetVector[i]->GetCubicVolume();
172 fNumTetMap[fMaterialVector[i]]++;
173 }
174 else
175 {
176 fVolumeMap[fMaterialVector[i]] = fTetVector[i]->GetCubicVolume();
177 fNumTetMap[fMaterialVector[i]] = 1;
178 }
179 }
180 ifpEle.close();
181 }
182
183 void TETModelImport::MaterialRead(G4String materialFile)
184 {
185 // Read material file (*.material)
186 //
187 std::ifstream ifpMat;
188
189 ifpMat.open((fPhantomDataPath + "/" + materialFile).c_str());
190 if(!ifpMat.is_open())
191 {
192 // exception for the case when there is no *.material file
193 G4Exception("TETModelImport::DataRead","",FatalErrorInArgument,
194 G4String(" There is no " + materialFile ).c_str());
195 }
196
197 G4cout << " Opening material file '" << materialFile << "'" <<G4endl;
198
199 char read_data[50];
200 char* token;
201 G4double zaid;
202 G4double fraction;
203 G4String MaterialName;
204 G4double density;
205 G4int i=0;
206
207 while(!ifpMat.eof())
208 {
209 ifpMat >> read_data; //ex) 'C' RBM
210 ifpMat >> MaterialName; //ex) C 'RBM'
211 ifpMat >> read_data;
212 density = std::atof(read_data); //ex) 1.30
213 ifpMat >> read_data; //ex) g/cm3
214 ifpMat >> read_data;
215 token = std::strtok(read_data,"m");
216 G4int matID = std::atoi(token); //ex) m'10'
217 fMaterialIndex.push_back(matID);
218 fOrganNameMap[matID]= MaterialName;
219 fDensityMap[matID] = density*g/cm3;
220
221 for(i=0 ; ; ++i)
222 {
223 ifpMat >> read_data;
224 if(std::strcmp(read_data, "C")==0 || ifpMat.eof()) break;
225
226 zaid = (G4int)(std::atoi(read_data)/1000);
227 ifpMat >> read_data;
228 fraction = -1.0 * std::atof(read_data);
229 fMaterialIndexMap[matID].push_back(std::make_pair(G4int(zaid), fraction));
230 }
231 }
232 ifpMat.close();
233
234 // Construct materials for each organ
235 //
236 G4NistManager* nistManager = G4NistManager::Instance();
237
238 for(i=0;i<(G4int)fMaterialIndex.size();++i)
239 {
240 G4int idx = fMaterialIndex[i];
241 auto* mat = new G4Material(fOrganNameMap[idx], fDensityMap[idx],
242 G4int(fMaterialIndexMap[idx].size()),
243 kStateSolid, NTP_Temperature, STP_Pressure);
244 for(G4int j=0;j<G4int(fMaterialIndexMap[idx].size());++j)
245 {
246 mat->AddElement(nistManager->FindOrBuildElement(fMaterialIndexMap[idx][j].first),
247 fMaterialIndexMap[idx][j].second);
248 }
249 fMaterialMap[idx]=mat;
250 fMassMap[idx]=fDensityMap[idx]*fVolumeMap[idx];
251 }
252 }
253
254 void TETModelImport::ColourRead()
255 {
256 // Read colour data file (colour.dat)
257 //
258 std::ifstream ifpColour;
259
260 ifpColour.open((fPhantomDataPath + "/" + "colour.dat").c_str());
261 if(!ifpColour.is_open())
262 {
263 // exception for the case when there is no colour.dat file
264 G4Exception("TETModelImport::DataRead","",FatalErrorInArgument,
265 G4String("Colour data file was not found ").c_str());
266 }
267
268 G4cout << " Opening colour data file 'colour.dat'" <<G4endl;
269
270 G4int organID;
271 G4double red, green, blue, alpha;
272 while( ifpColour >> organID >> red >> green >> blue >> alpha )
273 {
274 fColourMap[organID] = G4Colour(red, green, blue, alpha);
275 }
276
277 ifpColour.close();
278 }
279
280 void TETModelImport::PrintMaterialInfomation()
281 {
282 // Print the overall information for each organ
283 //
284 G4cout << G4endl
285 << std::setw(9) << "Organ ID"
286 << std::setw(11) << "# of Tet"
287 << std::setw(11) << "vol [cm3]"
288 << std::setw(11) << "d [g/cm3]"
289 << std::setw(11) << "mass [g]"
290 << "\t" << "organ/tissue"<< G4endl ;
291
292 G4cout << "--------------------------------------------------------------------------------"<<G4endl;
293
294 std::map<G4int, G4Material*>::const_iterator matIter;
295 G4cout<<std::setiosflags(std::ios::fixed);
296 G4cout.precision(3);
297 for(matIter=fMaterialMap.cbegin(); matIter!=fMaterialMap.cend(); ++matIter)
298 {
299 G4int idx = matIter->first;
300 G4cout << std::setw(9) << idx // organ ID
301 << std::setw(11) << fNumTetMap[idx] // # of tetrahedrons
302 << std::setw(11) << fVolumeMap[idx]/cm3 // organ volume
303 << std::setw(11) << fMaterialMap[idx] ->GetDensity()/(g/cm3) // organ density
304 << std::setw(11) << fMassMap[idx]/g // organ mass
305 << "\t"<<fMaterialMap[idx]->GetName() << G4endl; // organ name
306 }
307 }
308