| Geant4 Cross Reference |
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 // Author: Haegin Han
28 // Reference: ICRP Publication 145. Ann. ICRP
29 // Geant4 Contributors: J. Allison and S. Guat
30 //
31
32 #include "TETModelImport.hh"
33
34 TETModelImport::TETModelImport(G4bool isAF, G4
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_PAT
41 G4Exception("TETModelImport::TETModelImpor
42 G4String("PHANTOM_PATH environ
43 // default path for phantom data
44 fPhantomDataPath = "./ICRP145data";
45 }
46 else
47 {
48 // set path for phantom data as PHANTOM_PA
49 fPhantomDataPath = pPATH;
50 }
51
52 // set phantom name
53 if(!isAF) fPhantomName = "MRCP_AM";
54 else fPhantomName = "MRCP_AF";
55
56 G4cout << "=================================
57 G4cout << "\t" << fPhantomName << " was impl
58 G4cout << "=================================
59
60 G4String eleFile = fPhantomName + ".el
61 G4String nodeFile = fPhantomName + ".no
62 G4String materialFile = fPhantomName + ".ma
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 thi
69 if(ui) ColourRead();
70 // print the summary of phantom information
71 PrintMaterialInfomation();
72 }
73
74 void TETModelImport::DataRead(G4String eleFile
75 {
76 G4String tempStr;
77 G4int tempInt;
78
79 // Read *.node file
80 //
81 std::ifstream ifpNode;
82
83 ifpNode.open((fPhantomDataPath + "/" + nodeF
84 if(!ifpNode.is_open())
85 {
86 // exception for the case when there is no
87 G4Exception("TETModelImport::DataRead","",
88 G4String(" There is no "
89 }
90 G4cout << " Opening TETGEN node (vertex poi
91 << nodeFile << "'" <<G4endl;
92
93 G4int numVertex;
94 G4double xPos, yPos, zPos;
95 G4double xMin(DBL_MAX), yMin(DBL_MAX), zMin(
96 G4double xMax(DBL_MIN), yMax(DBL_MIN), zMax(
97
98 ifpNode >> numVertex >> 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::
110 fVertexVector.push_back(G4ThreeVector(xPos
111
112 // to get the information of the bounding
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 an
122 fBoundingBox_Min = G4ThreeVector(xMin,yMin,z
123 fBoundingBox_Max = G4ThreeVector(xMax,yMax,z
124 fPhantomSize = G4ThreeVector(xMax-xMin,yMax-
125
126 ifpNode.close();
127
128 // Read *.ele file
129 //
130 std::ifstream ifpEle;
131
132 ifpEle.open((fPhantomDataPath + "/" + eleFil
133 if(!ifpEle.is_open())
134 {
135 // exception for the case when there is no
136 G4Exception("TETModelImport::DataRead","",
137 G4String(" There is no "
138 }
139 G4cout << " Opening TETGEN elements (tetrah
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
159 fTetVector.push_back(new G4Tet("Tet_Solid"
160 fVertexVect
161 fVertexVect
162 fVertexVect
163 fVertexVect
164
165 // calculate the total volume and the numb
166 //std::map<G4int, G4double>::iterator Find
167 auto FindIter = fVolumeMap.find(fMaterialV
168
169 if(FindIter!=fVolumeMap.end())
170 {
171 FindIter->second += fTetVector[i]->GetCu
172 fNumTetMap[fMaterialVector[i]]++;
173 }
174 else
175 {
176 fVolumeMap[fMaterialVector[i]] = fTetVec
177 fNumTetMap[fMaterialVector[i]] = 1;
178 }
179 }
180 ifpEle.close();
181 }
182
183 void TETModelImport::MaterialRead(G4String mat
184 {
185 // Read material file (*.material)
186 //
187 std::ifstream ifpMat;
188
189 ifpMat.open((fPhantomDataPath + "/" + materi
190 if(!ifpMat.is_open())
191 {
192 // exception for the case when there is no
193 G4Exception("TETModelImport::DataRead","",
194 G4String(" There is no " + materialFi
195 }
196
197 G4cout << " Opening material file '" << mat
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; //e
210 ifpMat >> MaterialName; //e
211 ifpMat >> read_data;
212 density = std::atof(read_data); //e
213 ifpMat >> read_data; //e
214 ifpMat >> read_data;
215 token = std::strtok(read_data,"m");
216 G4int matID = std::atoi(token); //e
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 || ifp
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::
230 }
231 }
232 ifpMat.close();
233
234 // Construct materials for each organ
235 //
236 G4NistManager* nistManager = G4NistManager::
237
238 for(i=0;i<(G4int)fMaterialIndex.size();++i)
239 {
240 G4int idx = fMaterialIndex[i];
241 auto* mat = new G4Material(fOrganNameMap[i
242 G4int(fMaterial
243 kStateSolid, NT
244 for(G4int j=0;j<G4int(fMaterialIndexMap[id
245 {
246 mat->AddElement(nistManager->FindOrBuild
247
248 }
249 fMaterialMap[idx]=mat;
250 fMassMap[idx]=fDensityMap[idx]*fVolumeMap[
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 + "/" + "co
261 if(!ifpColour.is_open())
262 {
263 // exception for the case when there is no
264 G4Exception("TETModelImport::DataRead","",
265 G4String("Colour data file was
266 }
267
268 G4cout << " Opening colour data file 'colou
269
270 G4int organID;
271 G4double red, green, blue, alpha;
272 while( ifpColour >> organID >> red >> green
273 {
274 fColourMap[organID] = G4Colour(red, green,
275 }
276
277 ifpColour.close();
278 }
279
280 void TETModelImport::PrintMaterialInfomation()
281 {
282 // Print the overall information for each or
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 << "---------------------------------
293
294 std::map<G4int, G4Material*>::const_iterator
295 G4cout<<std::setiosflags(std::ios::fixed);
296 G4cout.precision(3);
297 for(matIter=fMaterialMap.cbegin(); matIter!=
298 {
299 G4int idx = matIter->first;
300 G4cout << std::setw(9) << idx
301 << std::setw(11) << fNumTetMap[idx]
302 << std::setw(11) << fVolumeMap[idx]
303 << std::setw(11) << fMaterialMap[id
304 << std::setw(11) << fMassMap[idx]/g
305 << "\t"<<fMaterialMap[idx]->GetName
306 }
307 }
308