| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 // Author: Haegin Han 27 // Author: Haegin Han
28 // Reference: ICRP Publication 145. Ann. ICRP 28 // Reference: ICRP Publication 145. Ann. ICRP 49(3), 2020.
29 // Geant4 Contributors: J. Allison and S. Guat 29 // Geant4 Contributors: J. Allison and S. Guatelli
30 // 30 //
31 31
32 #include "TETModelImport.hh" 32 #include "TETModelImport.hh"
33 33
34 TETModelImport::TETModelImport(G4bool isAF, G4 34 TETModelImport::TETModelImport(G4bool isAF, G4UIExecutive* ui)
35 { 35 {
36 // set path for phantom data 36 // set path for phantom data
37 char* pPATH = std::getenv("PHANTOM_PATH"); 37 char* pPATH = std::getenv("PHANTOM_PATH");
38 if( pPATH == nullptr ) 38 if( pPATH == nullptr )
39 { 39 {
40 // exception for the case when PHANTOM_PAT 40 // exception for the case when PHANTOM_PATH environment variable was not set
41 G4Exception("TETModelImport::TETModelImpor 41 G4Exception("TETModelImport::TETModelImport","",JustWarning,
42 G4String("PHANTOM_PATH environ 42 G4String("PHANTOM_PATH environment variable was not set. Using ./ICRP145data.").c_str());
43 // default path for phantom data 43 // default path for phantom data
44 fPhantomDataPath = "./ICRP145data"; 44 fPhantomDataPath = "./ICRP145data";
45 } 45 }
46 else 46 else
47 { 47 {
48 // set path for phantom data as PHANTOM_PA 48 // set path for phantom data as PHANTOM_PATH
49 fPhantomDataPath = pPATH; 49 fPhantomDataPath = pPATH;
50 } 50 }
51 51
52 // set phantom name 52 // set phantom name
53 if(!isAF) fPhantomName = "MRCP_AM"; 53 if(!isAF) fPhantomName = "MRCP_AM";
54 else fPhantomName = "MRCP_AF"; 54 else fPhantomName = "MRCP_AF";
55 55
56 G4cout << "================================= 56 G4cout << "================================================================================"<<G4endl;
57 G4cout << "\t" << fPhantomName << " was impl 57 G4cout << "\t" << fPhantomName << " was implemented in this CODE!! "<< G4endl;
58 G4cout << "================================= 58 G4cout << "================================================================================"<<G4endl;
59 59
60 G4String eleFile = fPhantomName + ".el 60 G4String eleFile = fPhantomName + ".ele";
61 G4String nodeFile = fPhantomName + ".no 61 G4String nodeFile = fPhantomName + ".node";
62 G4String materialFile = fPhantomName + ".ma 62 G4String materialFile = fPhantomName + ".material";
63 63
64 // read phantom data files (*. ele, *.node) 64 // read phantom data files (*. ele, *.node)
65 DataRead(eleFile, nodeFile); 65 DataRead(eleFile, nodeFile);
66 // read material file (*.material) 66 // read material file (*.material)
67 MaterialRead(materialFile); 67 MaterialRead(materialFile);
68 // read colour data file (colour.dat) if thi 68 // read colour data file (colour.dat) if this is interactive mode
69 if(ui) ColourRead(); 69 if(ui) ColourRead();
70 // print the summary of phantom information 70 // print the summary of phantom information
71 PrintMaterialInfomation(); 71 PrintMaterialInfomation();
72 } 72 }
73 73
74 void TETModelImport::DataRead(G4String eleFile 74 void TETModelImport::DataRead(G4String eleFile, G4String nodeFile)
75 { 75 {
76 G4String tempStr; 76 G4String tempStr;
77 G4int tempInt; 77 G4int tempInt;
78 78
79 // Read *.node file 79 // Read *.node file
80 // 80 //
81 std::ifstream ifpNode; 81 std::ifstream ifpNode;
82 82
83 ifpNode.open((fPhantomDataPath + "/" + nodeF 83 ifpNode.open((fPhantomDataPath + "/" + nodeFile).c_str());
84 if(!ifpNode.is_open()) 84 if(!ifpNode.is_open())
85 { 85 {
86 // exception for the case when there is no 86 // exception for the case when there is no *.node file
87 G4Exception("TETModelImport::DataRead","", 87 G4Exception("TETModelImport::DataRead","",FatalErrorInArgument,
88 G4String(" There is no " 88 G4String(" There is no " + nodeFile ).c_str());
89 } 89 }
90 G4cout << " Opening TETGEN node (vertex poi 90 G4cout << " Opening TETGEN node (vertex points: x y z) file '"
91 << nodeFile << "'" <<G4endl; 91 << nodeFile << "'" <<G4endl;
92 92
93 G4int numVertex; 93 G4int numVertex;
94 G4double xPos, yPos, zPos; 94 G4double xPos, yPos, zPos;
95 G4double xMin(DBL_MAX), yMin(DBL_MAX), zMin( 95 G4double xMin(DBL_MAX), yMin(DBL_MAX), zMin(DBL_MAX);
96 G4double xMax(DBL_MIN), yMax(DBL_MIN), zMax( 96 G4double xMax(DBL_MIN), yMax(DBL_MIN), zMax(DBL_MIN);
97 97
98 ifpNode >> numVertex >> tempInt >> tempInt > 98 ifpNode >> numVertex >> tempInt >> tempInt >> tempInt;
99 99
100 for(G4int i=0; i<numVertex; ++i) 100 for(G4int i=0; i<numVertex; ++i)
101 { 101 {
102 ifpNode >> tempInt >> xPos >> yPos >> zPos 102 ifpNode >> tempInt >> xPos >> yPos >> zPos;
103 103
104 // set the unit 104 // set the unit
105 xPos*=cm; 105 xPos*=cm;
106 yPos*=cm; 106 yPos*=cm;
107 zPos*=cm; 107 zPos*=cm;
108 108
109 // save the node data as the form of std:: 109 // save the node data as the form of std::vector<G4ThreeVector>
110 fVertexVector.push_back(G4ThreeVector(xPos 110 fVertexVector.push_back(G4ThreeVector(xPos, yPos, zPos));
111 111
112 // to get the information of the bounding 112 // to get the information of the bounding box of phantom
113 if (xPos < xMin) xMin = xPos; 113 if (xPos < xMin) xMin = xPos;
114 if (xPos > xMax) xMax = xPos; 114 if (xPos > xMax) xMax = xPos;
115 if (yPos < yMin) yMin = yPos; 115 if (yPos < yMin) yMin = yPos;
116 if (yPos > yMax) yMax = yPos; 116 if (yPos > yMax) yMax = yPos;
117 if (zPos < zMin) zMin = zPos; 117 if (zPos < zMin) zMin = zPos;
118 if (zPos > zMax) zMax = zPos; 118 if (zPos > zMax) zMax = zPos;
119 } 119 }
120 120
121 // set the variables for the bounding box an 121 // set the variables for the bounding box and phantom size
122 fBoundingBox_Min = G4ThreeVector(xMin,yMin,z 122 fBoundingBox_Min = G4ThreeVector(xMin,yMin,zMin);
123 fBoundingBox_Max = G4ThreeVector(xMax,yMax,z 123 fBoundingBox_Max = G4ThreeVector(xMax,yMax,zMax);
124 fPhantomSize = G4ThreeVector(xMax-xMin,yMax- 124 fPhantomSize = G4ThreeVector(xMax-xMin,yMax-yMin,zMax-zMin);
125 125
126 ifpNode.close(); 126 ifpNode.close();
127 127
128 // Read *.ele file 128 // Read *.ele file
129 // 129 //
130 std::ifstream ifpEle; 130 std::ifstream ifpEle;
131 131
132 ifpEle.open((fPhantomDataPath + "/" + eleFil 132 ifpEle.open((fPhantomDataPath + "/" + eleFile).c_str());
133 if(!ifpEle.is_open()) 133 if(!ifpEle.is_open())
134 { 134 {
135 // exception for the case when there is no 135 // exception for the case when there is no *.ele file
136 G4Exception("TETModelImport::DataRead","", 136 G4Exception("TETModelImport::DataRead","",FatalErrorInArgument,
137 G4String(" There is no " 137 G4String(" There is no " + eleFile ).c_str());
138 } 138 }
139 G4cout << " Opening TETGEN elements (tetrah 139 G4cout << " Opening TETGEN elements (tetrahedron with node No.) file '"
140 << eleFile << "'" <<G4endl; 140 << eleFile << "'" <<G4endl;
141 141
142 G4int numEle; 142 G4int numEle;
143 ifpEle >> numEle >> tempInt >> tempInt; 143 ifpEle >> numEle >> tempInt >> tempInt;
144 144
145 for(G4int i=0; i<numEle; ++i) 145 for(G4int i=0; i<numEle; ++i)
146 { 146 {
147 ifpEle >> tempInt; 147 ifpEle >> tempInt;
148 auto* ele = new G4int[4]; 148 auto* ele = new G4int[4];
149 for(G4int j=0;j<4;j++) 149 for(G4int j=0;j<4;j++)
150 { 150 {
151 ifpEle >> tempInt; 151 ifpEle >> tempInt;
152 ele[j]=tempInt; 152 ele[j]=tempInt;
153 } 153 }
154 fEleVector.push_back(ele); 154 fEleVector.push_back(ele);
155 ifpEle >> tempInt; 155 ifpEle >> tempInt;
156 fMaterialVector.push_back(tempInt); 156 fMaterialVector.push_back(tempInt);
157 157
158 // save the element (tetrahedron) data as 158 // save the element (tetrahedron) data as the form of std::vector<G4Tet*>
159 fTetVector.push_back(new G4Tet("Tet_Solid" 159 fTetVector.push_back(new G4Tet("Tet_Solid",
160 fVertexVect 160 fVertexVector[ele[0]],
161 fVertexVect 161 fVertexVector[ele[1]],
162 fVertexVect 162 fVertexVector[ele[2]],
163 fVertexVect 163 fVertexVector[ele[3]]));
164 164
165 // calculate the total volume and the numb 165 // calculate the total volume and the number of tetrahedrons for each organ
166 //std::map<G4int, G4double>::iterator Find 166 //std::map<G4int, G4double>::iterator FindIter = fVolumeMap.find(fMaterialVector[i]);
167 auto FindIter = fVolumeMap.find(fMaterialV 167 auto FindIter = fVolumeMap.find(fMaterialVector[i]);
168 168
169 if(FindIter!=fVolumeMap.end()) 169 if(FindIter!=fVolumeMap.end())
170 { 170 {
171 FindIter->second += fTetVector[i]->GetCu 171 FindIter->second += fTetVector[i]->GetCubicVolume();
172 fNumTetMap[fMaterialVector[i]]++; 172 fNumTetMap[fMaterialVector[i]]++;
173 } 173 }
174 else 174 else
175 { 175 {
176 fVolumeMap[fMaterialVector[i]] = fTetVec 176 fVolumeMap[fMaterialVector[i]] = fTetVector[i]->GetCubicVolume();
177 fNumTetMap[fMaterialVector[i]] = 1; 177 fNumTetMap[fMaterialVector[i]] = 1;
178 } 178 }
179 } 179 }
180 ifpEle.close(); 180 ifpEle.close();
181 } 181 }
182 182
183 void TETModelImport::MaterialRead(G4String mat 183 void TETModelImport::MaterialRead(G4String materialFile)
184 { 184 {
185 // Read material file (*.material) 185 // Read material file (*.material)
186 // 186 //
187 std::ifstream ifpMat; 187 std::ifstream ifpMat;
188 188
189 ifpMat.open((fPhantomDataPath + "/" + materi 189 ifpMat.open((fPhantomDataPath + "/" + materialFile).c_str());
190 if(!ifpMat.is_open()) 190 if(!ifpMat.is_open())
191 { 191 {
192 // exception for the case when there is no 192 // exception for the case when there is no *.material file
193 G4Exception("TETModelImport::DataRead","", 193 G4Exception("TETModelImport::DataRead","",FatalErrorInArgument,
194 G4String(" There is no " + materialFi 194 G4String(" There is no " + materialFile ).c_str());
195 } 195 }
196 196
197 G4cout << " Opening material file '" << mat 197 G4cout << " Opening material file '" << materialFile << "'" <<G4endl;
198 198
199 char read_data[50]; 199 char read_data[50];
200 char* token; 200 char* token;
201 G4double zaid; 201 G4double zaid;
202 G4double fraction; 202 G4double fraction;
203 G4String MaterialName; 203 G4String MaterialName;
204 G4double density; 204 G4double density;
205 G4int i=0; 205 G4int i=0;
206 206
207 while(!ifpMat.eof()) 207 while(!ifpMat.eof())
208 { 208 {
209 ifpMat >> read_data; //e 209 ifpMat >> read_data; //ex) 'C' RBM
210 ifpMat >> MaterialName; //e 210 ifpMat >> MaterialName; //ex) C 'RBM'
211 ifpMat >> read_data; 211 ifpMat >> read_data;
212 density = std::atof(read_data); //e 212 density = std::atof(read_data); //ex) 1.30
213 ifpMat >> read_data; //e 213 ifpMat >> read_data; //ex) g/cm3
214 ifpMat >> read_data; 214 ifpMat >> read_data;
215 token = std::strtok(read_data,"m"); 215 token = std::strtok(read_data,"m");
216 G4int matID = std::atoi(token); //e 216 G4int matID = std::atoi(token); //ex) m'10'
217 fMaterialIndex.push_back(matID); 217 fMaterialIndex.push_back(matID);
218 fOrganNameMap[matID]= MaterialName; 218 fOrganNameMap[matID]= MaterialName;
219 fDensityMap[matID] = density*g/cm3; 219 fDensityMap[matID] = density*g/cm3;
220 220
221 for(i=0 ; ; ++i) 221 for(i=0 ; ; ++i)
222 { 222 {
223 ifpMat >> read_data; 223 ifpMat >> read_data;
224 if(std::strcmp(read_data, "C")==0 || ifp 224 if(std::strcmp(read_data, "C")==0 || ifpMat.eof()) break;
225 225
226 zaid = (G4int)(std::atoi(read_data)/1000 226 zaid = (G4int)(std::atoi(read_data)/1000);
227 ifpMat >> read_data; 227 ifpMat >> read_data;
228 fraction = -1.0 * std::atof(read_data); 228 fraction = -1.0 * std::atof(read_data);
229 fMaterialIndexMap[matID].push_back(std:: 229 fMaterialIndexMap[matID].push_back(std::make_pair(G4int(zaid), fraction));
230 } 230 }
231 } 231 }
232 ifpMat.close(); 232 ifpMat.close();
233 233
234 // Construct materials for each organ 234 // Construct materials for each organ
235 // 235 //
236 G4NistManager* nistManager = G4NistManager:: 236 G4NistManager* nistManager = G4NistManager::Instance();
237 237
238 for(i=0;i<(G4int)fMaterialIndex.size();++i) 238 for(i=0;i<(G4int)fMaterialIndex.size();++i)
239 { 239 {
240 G4int idx = fMaterialIndex[i]; 240 G4int idx = fMaterialIndex[i];
241 auto* mat = new G4Material(fOrganNameMap[i 241 auto* mat = new G4Material(fOrganNameMap[idx], fDensityMap[idx],
242 G4int(fMaterial 242 G4int(fMaterialIndexMap[idx].size()),
243 kStateSolid, NT 243 kStateSolid, NTP_Temperature, STP_Pressure);
244 for(G4int j=0;j<G4int(fMaterialIndexMap[id 244 for(G4int j=0;j<G4int(fMaterialIndexMap[idx].size());++j)
245 { 245 {
246 mat->AddElement(nistManager->FindOrBuild 246 mat->AddElement(nistManager->FindOrBuildElement(fMaterialIndexMap[idx][j].first),
247 247 fMaterialIndexMap[idx][j].second);
248 } 248 }
249 fMaterialMap[idx]=mat; 249 fMaterialMap[idx]=mat;
250 fMassMap[idx]=fDensityMap[idx]*fVolumeMap[ 250 fMassMap[idx]=fDensityMap[idx]*fVolumeMap[idx];
251 } 251 }
252 } 252 }
253 253
254 void TETModelImport::ColourRead() 254 void TETModelImport::ColourRead()
255 { 255 {
256 // Read colour data file (colour.dat) 256 // Read colour data file (colour.dat)
257 // 257 //
258 std::ifstream ifpColour; 258 std::ifstream ifpColour;
259 259
260 ifpColour.open((fPhantomDataPath + "/" + "co 260 ifpColour.open((fPhantomDataPath + "/" + "colour.dat").c_str());
261 if(!ifpColour.is_open()) 261 if(!ifpColour.is_open())
262 { 262 {
263 // exception for the case when there is no 263 // exception for the case when there is no colour.dat file
264 G4Exception("TETModelImport::DataRead","", 264 G4Exception("TETModelImport::DataRead","",FatalErrorInArgument,
265 G4String("Colour data file was 265 G4String("Colour data file was not found ").c_str());
266 } 266 }
267 267
268 G4cout << " Opening colour data file 'colou 268 G4cout << " Opening colour data file 'colour.dat'" <<G4endl;
269 269
270 G4int organID; 270 G4int organID;
271 G4double red, green, blue, alpha; 271 G4double red, green, blue, alpha;
272 while( ifpColour >> organID >> red >> green 272 while( ifpColour >> organID >> red >> green >> blue >> alpha )
273 { 273 {
274 fColourMap[organID] = G4Colour(red, green, 274 fColourMap[organID] = G4Colour(red, green, blue, alpha);
275 } 275 }
276 276
277 ifpColour.close(); 277 ifpColour.close();
278 } 278 }
279 279
280 void TETModelImport::PrintMaterialInfomation() 280 void TETModelImport::PrintMaterialInfomation()
281 { 281 {
282 // Print the overall information for each or 282 // Print the overall information for each organ
283 // 283 //
284 G4cout << G4endl 284 G4cout << G4endl
285 << std::setw(9) << "Organ ID" 285 << std::setw(9) << "Organ ID"
286 << std::setw(11) << "# of Tet" 286 << std::setw(11) << "# of Tet"
287 << std::setw(11) << "vol [cm3]" 287 << std::setw(11) << "vol [cm3]"
288 << std::setw(11) << "d [g/cm3]" 288 << std::setw(11) << "d [g/cm3]"
289 << std::setw(11) << "mass [g]" 289 << std::setw(11) << "mass [g]"
290 << "\t" << "organ/tissue"<< G4endl ; 290 << "\t" << "organ/tissue"<< G4endl ;
291 291
292 G4cout << "--------------------------------- 292 G4cout << "--------------------------------------------------------------------------------"<<G4endl;
293 293
294 std::map<G4int, G4Material*>::const_iterator 294 std::map<G4int, G4Material*>::const_iterator matIter;
295 G4cout<<std::setiosflags(std::ios::fixed); 295 G4cout<<std::setiosflags(std::ios::fixed);
296 G4cout.precision(3); 296 G4cout.precision(3);
297 for(matIter=fMaterialMap.cbegin(); matIter!= 297 for(matIter=fMaterialMap.cbegin(); matIter!=fMaterialMap.cend(); ++matIter)
298 { 298 {
299 G4int idx = matIter->first; 299 G4int idx = matIter->first;
300 G4cout << std::setw(9) << idx 300 G4cout << std::setw(9) << idx // organ ID
301 << std::setw(11) << fNumTetMap[idx] 301 << std::setw(11) << fNumTetMap[idx] // # of tetrahedrons
302 << std::setw(11) << fVolumeMap[idx] 302 << std::setw(11) << fVolumeMap[idx]/cm3 // organ volume
303 << std::setw(11) << fMaterialMap[id 303 << std::setw(11) << fMaterialMap[idx] ->GetDensity()/(g/cm3) // organ density
304 << std::setw(11) << fMassMap[idx]/g 304 << std::setw(11) << fMassMap[idx]/g // organ mass
305 << "\t"<<fMaterialMap[idx]->GetName 305 << "\t"<<fMaterialMap[idx]->GetName() << G4endl; // organ name
306 } 306 }
307 } 307 }
308 308