| 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 // <<
27 /// \file DicomPhantomZSliceHeader.cc 26 /// \file DicomPhantomZSliceHeader.cc
28 /// \brief Implementation of the DicomPhantomZ 27 /// \brief Implementation of the DicomPhantomZSliceHeader class
29 // <<
30 28
31 #include "DicomPhantomZSliceHeader.hh" << 29 #include "globals.hh"
32 <<
33 #include "G4GeometryTolerance.hh" <<
34 #include "G4LogicalVolume.hh" 30 #include "G4LogicalVolume.hh"
35 #include "G4Material.hh" <<
36 #include "G4MaterialTable.hh" 31 #include "G4MaterialTable.hh"
37 #include "G4NistManager.hh" << 32 #include "G4Material.hh"
38 #include "globals.hh" << 33 #include "G4GeometryTolerance.hh"
39 34
40 //....oooOO0OOooo........oooOO0OOooo........oo << 35 #include "DicomPhantomZSliceHeader.hh"
41 DicomPhantomZSliceHeader::DicomPhantomZSliceHe <<
42 : fNoVoxelsX(0), <<
43 fNoVoxelsY(0), <<
44 fNoVoxelsZ(0), <<
45 fMinX(0), <<
46 fMinY(0), <<
47 fMinZ(0), <<
48 fMaxX(0), <<
49 fMaxY(0), <<
50 fMaxZ(0), <<
51 fFilename(fname), <<
52 fSliceLocation(0) <<
53 {} <<
54 36
55 //....oooOO0OOooo........oooOO0OOooo........oo << 37 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56 DicomPhantomZSliceHeader::~DicomPhantomZSliceH << 38 DicomPhantomZSliceHeader::DicomPhantomZSliceHeader( const DicomPhantomZSliceHeader& rhs )
>> 39 {
>> 40 fNoVoxelX = rhs.GetNoVoxelX();
>> 41 fNoVoxelY = rhs.GetNoVoxelY();
>> 42 fNoVoxelZ = rhs.GetNoVoxelZ();
>> 43 fMinX = rhs.GetMinX();
>> 44 fMaxX = rhs.GetMaxX();
>> 45 fMinY = rhs.GetMinY();
>> 46 fMaxY = rhs.GetMaxY();
>> 47 fMinZ = rhs.GetMinZ();
>> 48 fMaxZ = rhs.GetMaxZ();
>> 49 fMaterialNames = rhs.GetMaterialNames();
57 50
58 //....oooOO0OOooo........oooOO0OOooo........oo << 51 }
59 DicomPhantomZSliceHeader::DicomPhantomZSliceHe << 52
>> 53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 54 DicomPhantomZSliceHeader::DicomPhantomZSliceHeader( std::ifstream& fin )
60 { 55 {
61 //----- Read material indices and names 56 //----- Read material indices and names
62 G4int nmate; 57 G4int nmate;
63 G4String mateindex; 58 G4String mateindex;
64 G4String matename; 59 G4String matename;
65 fin >> nmate; 60 fin >> nmate;
66 #ifdef G4VERBOSE 61 #ifdef G4VERBOSE
67 G4cout << " DicomPhantomZSliceHeader reading 62 G4cout << " DicomPhantomZSliceHeader reading number of materials " << nmate << G4endl;
68 #endif 63 #endif
69 64
70 for (G4int im = 0; im < nmate; im++) { << 65 for( G4int im = 0; im < nmate; im++ ){
71 fin >> mateindex >> matename; 66 fin >> mateindex >> matename;
72 #ifdef G4VERBOSE 67 #ifdef G4VERBOSE
73 // G4cout << " DicomPhantomZSliceHeader re << 68 G4cout << " DicomPhantomZSliceHeader reading material " << im << " : " << mateindex << " " << matename << G4endl;
74 // << im << " : "<< mateindex << " " << <<
75 #endif 69 #endif
76 70
77 if (!CheckMaterialExists(matename)) { << 71 if( ! CheckMaterialExists( matename ) ) {
78 G4Exception("DicomPhantomZSliceHeader::D << 72 G4Exception("DicomPhantomZSliceHeader::DicomPhantomZSliceHeader","A material is found in file that is not built in the C++ code",FatalErrorInArgument,matename.c_str());
79 "A material is found in file <<
80 FatalErrorInArgument, matena <<
81 } 73 }
82 74
83 fMaterialNames.push_back(matename); 75 fMaterialNames.push_back(matename);
84 } 76 }
85 77
86 //----- Read number of voxels 78 //----- Read number of voxels
87 fin >> fNoVoxelsX >> fNoVoxelsY >> fNoVoxels << 79 fin >> fNoVoxelX >> fNoVoxelY >> fNoVoxelZ;
88 #ifdef G4VERBOSE 80 #ifdef G4VERBOSE
89 G4cout << " Number of voxels " << fNoVoxelsX << 81 G4cout << " Number of voxels " << fNoVoxelX << " " << fNoVoxelY << " " << fNoVoxelZ << G4endl;
90 #endif 82 #endif
91 83
92 //----- Read minimal and maximal extensions 84 //----- Read minimal and maximal extensions (= walls of phantom)
93 fin >> fMinX >> fMaxX; 85 fin >> fMinX >> fMaxX;
94 fin >> fMinY >> fMaxY; 86 fin >> fMinY >> fMaxY;
95 fin >> fMinZ >> fMaxZ; 87 fin >> fMinZ >> fMaxZ;
96 #ifdef G4VERBOSE 88 #ifdef G4VERBOSE
97 G4cout << " Extension in X " << fMinX << " " << 89 G4cout << " Extension in X " << fMinX << " " << fMaxX << G4endl
98 << " " << fMaxY << G4endl << " Extens << 90 << " Extension in Y " << fMinY << " " << fMaxY << G4endl
>> 91 << " Extension in Z " << fMinZ << " " << fMaxZ << G4endl;
99 #endif 92 #endif
100 93
101 fSliceLocation = 0.5 * (fMinZ + fMaxZ); <<
102 } 94 }
103 95
104 //....oooOO0OOooo........oooOO0OOooo........oo << 96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
105 G4bool DicomPhantomZSliceHeader::CheckMaterial << 97 G4bool DicomPhantomZSliceHeader::CheckMaterialExists( const G4String& mateName )
106 { 98 {
>> 99 G4bool bFound = FALSE;
>> 100
107 const G4MaterialTable* matTab = G4Material:: 101 const G4MaterialTable* matTab = G4Material::GetMaterialTable();
108 std::vector<G4Material*>::const_iterator mat 102 std::vector<G4Material*>::const_iterator matite;
109 for (matite = matTab->begin(); matite != mat << 103 for( matite = matTab->begin(); matite != matTab->end(); matite++ ) {
110 if ((*matite)->GetName() == mateName) { << 104 if( (*matite)->GetName() == mateName ) {
111 return true; << 105 bFound = TRUE;
>> 106 break;
112 } 107 }
113 } 108 }
114 109
115 G4Material* g4mate = G4NistManager::Instance << 110 return bFound;
116 if (g4mate) { << 111
117 return false; <<
118 } <<
119 else { <<
120 return true; <<
121 } <<
122 } 112 }
123 113
124 //....oooOO0OOooo........oooOO0OOooo........oo << 114 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
125 void DicomPhantomZSliceHeader::operator+=(cons << 115 void DicomPhantomZSliceHeader::operator+=( const DicomPhantomZSliceHeader& rhs )
126 { 116 {
127 *this = *this + rhs; 117 *this = *this + rhs;
128 } 118 }
129 119
130 //....oooOO0OOooo........oooOO0OOooo........oo << 120 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
131 DicomPhantomZSliceHeader DicomPhantomZSliceHea << 121 DicomPhantomZSliceHeader DicomPhantomZSliceHeader::operator+( const DicomPhantomZSliceHeader& rhs )
132 { 122 {
133 //----- Check that both slices has the same 123 //----- Check that both slices has the same dimensions
134 if (fNoVoxelsX != rhs.GetNoVoxelsX() || fNoV << 124 if( fNoVoxelX != rhs.GetNoVoxelX()
135 G4cerr << "DicomPhantomZSliceHeader error << 125 || fNoVoxelY != rhs.GetNoVoxelY() ) {
136 !!! Different number of voxels: " << 126 G4cerr << "DicomPhantomZSliceHeader error adding two slice headers: !!! Different number of voxels: "
137 << " X= " << fNoVoxelsX << " =? " << 127 << " X= " << fNoVoxelX << " =? " << rhs.GetNoVoxelX()
138 << " =? " << rhs.GetNoVoxelsY() << << 128 << " Y= " << fNoVoxelY << " =? " << rhs.GetNoVoxelY()
>> 129 << " Z= " << fNoVoxelZ << " =? " << rhs.GetNoVoxelZ()
139 << G4endl; 130 << G4endl;
140 G4Exception("DicomPhantomZSliceHeader::Dic << 131 G4Exception("DicomPhantomZSliceHeader::DicomPhantomZSliceHeader","",FatalErrorInArgument,"");
141 } 132 }
142 //----- Check that both slices has the same 133 //----- Check that both slices has the same extensions
143 if (fMinX != rhs.GetMinX() || fMaxX != rhs.G << 134 if( fMinX != rhs.GetMinX() || fMaxX != rhs.GetMaxX()
144 || fMaxY != rhs.GetMaxY()) << 135 || fMinY != rhs.GetMinY() || fMaxY != rhs.GetMaxY() ) {
145 { << 136 G4cerr << "DicomPhantomZSliceHeader error adding two slice headers: !!! Different extensions: "
146 G4cerr << "DicomPhantomZSliceHeader error << 137 << " Xmin= " << fMinX << " =? " << rhs.GetMinX()
147 !!! Different extensions: " << 138 << " Xmax= " << fMaxX << " =? " << rhs.GetMaxX()
148 << " Xmin= " << fMinX << " =? " << << 139 << " Ymin= " << fMinY << " =? " << rhs.GetMinY()
149 << rhs.GetMaxX() << " Ymin= " << f << 140 << " Ymax= " << fMaxY << " =? " << rhs.GetMaxY()
150 << " =? " << rhs.GetMaxY() << G4end << 141 << G4endl;
151 G4Exception("DicomPhantomZSliceHeader::ope << 142 G4Exception("DicomPhantomZSliceHeader::operator+","",FatalErrorInArgument,"");
152 } 143 }
153 << 144
154 //----- Check that both slices have the same << 145 //----- Check that both slices has the same materials
155 std::vector<G4String> fMaterialNames2 = rhs. 146 std::vector<G4String> fMaterialNames2 = rhs.GetMaterialNames();
156 if (fMaterialNames.size() != fMaterialNames2 << 147 if( fMaterialNames.size() != fMaterialNames2.size() ) {
157 G4cerr << "DicomPhantomZSliceHeader error << 148 G4cerr << "DicomPhantomZSliceHeader error adding two slice headers: !!! Different number of materials: " << fMaterialNames.size() << " =? " << fMaterialNames2.size() << G4endl;
158 !!! Different number of materials: " << 149 G4Exception("DicomPhantomZSliceHeader::operator+","",FatalErrorInArgument,"");
159 << fMaterialNames.size() << " =? " <<
160 G4Exception("DicomPhantomZSliceHeader::ope <<
161 } 150 }
162 for (unsigned int ii = 0; ii < fMaterialName << 151 for( unsigned int ii = 0; ii < fMaterialNames.size(); ii++ ) {
163 if (fMaterialNames[ii] != fMaterialNames2[ << 152 if( fMaterialNames[ii] != fMaterialNames2[ii] ) {
164 G4cerr << "DicomPhantomZSliceHeader erro << 153 G4cerr << "DicomPhantomZSliceHeader error adding two slice headers: !!! Different material number " << ii << " : " << fMaterialNames[ii] << " =? " << fMaterialNames2[ii] << G4endl;
165 !!! Different material number " << 154 G4Exception("DicomPhantomZSliceHeader::operator+","",FatalErrorInArgument,"");
166 << ii << " : " << fMaterialNames[ <<
167 G4Exception("DicomPhantomZSliceHeader::o <<
168 } 155 }
169 } 156 }
170 << 157
171 //----- Check that the slices are contiguous 158 //----- Check that the slices are contiguous in Z
172 if (std::fabs(fMinZ - rhs.GetMaxZ()) > G4Geo << 159 if( std::fabs( fMinZ - rhs.GetMaxZ() ) > G4GeometryTolerance::GetInstance()->GetRadialTolerance() &&
173 && std::fabs(fMaxZ - rhs.GetMinZ()) << 160 std::fabs( fMaxZ - rhs.GetMinZ() ) > G4GeometryTolerance::GetInstance()->GetRadialTolerance() ){
174 > G4GeometryTolerance::GetInstance( << 161 G4cerr << "DicomPhantomZSliceHeader error adding two slice headers: !!! Slices are not contiguous in Z "
175 { << 162 << " Zmin= " << fMinZ << " & " << rhs.GetMinZ()
176 G4cerr << "DicomPhantomZSliceHeader error << 163 << " Zmax= " << fMaxZ << " & " << rhs.GetMaxZ()
177 Slices are not contiguous in Z " << 164 << G4endl;
178 << " Zmin= " << fMinZ << " & " << << 165 G4Exception("DicomPhantomZSliceHeader::operator+","",FatalErrorInArgument,"");
179 << rhs.GetMaxZ() << G4endl; <<
180 G4Exception("DicomPhantomZSliceHeader::ope <<
181 } 166 }
182 167
183 //----- Build slice header copying first one 168 //----- Build slice header copying first one
184 DicomPhantomZSliceHeader temp(*this); << 169 DicomPhantomZSliceHeader temp( *this );
185 170
186 //----- Add data from second slice header 171 //----- Add data from second slice header
187 temp.SetMinZ(std::min(fMinZ, rhs.GetMinZ())) << 172 temp.SetMinZ( std::min( fMinZ, rhs.GetMinZ() ) );
188 temp.SetMaxZ(std::max(fMaxZ, rhs.GetMaxZ())) << 173 temp.SetMaxZ( std::max( fMaxZ, rhs.GetMaxZ() ) );
189 temp.SetNoVoxelsZ(fNoVoxelsZ + rhs.GetNoVoxe << 174 temp.SetNoVoxelZ( fNoVoxelZ + rhs.GetNoVoxelZ() );
190 175
191 return temp; 176 return temp;
192 } 177 }
193 <<
194 //....oooOO0OOooo........oooOO0OOooo........oo <<
195 void DicomPhantomZSliceHeader::DumpToFile() <<
196 { <<
197 G4cout << "DicomPhantomZSliceHeader::Dumping <<
198 // sleep(5); <<
199 <<
200 // May seen counter-intuitive (dumping to f <<
201 // the reason for this is modification slic <<
202 if (fMateIDs.size() == 0 || fValues.size() = <<
203 ReadDataFromFile(); <<
204 } <<
205 <<
206 std::ofstream out; <<
207 out.open(fFilename.c_str()); <<
208 <<
209 if (!out) { <<
210 G4String descript = "DicomPhantomZSliceHea <<
211 G4Exception(descript.c_str(), "", FatalExc <<
212 } <<
213 <<
214 out << fMaterialNames.size() << std::endl; <<
215 for (unsigned int i = 0; i < fMaterialNames. <<
216 out << i << " " << fMaterialNames.at(i) << <<
217 } <<
218 <<
219 out << fNoVoxelsX << " " << fNoVoxelsY << " <<
220 out << fMinX << " " << fMaxX << std::endl; <<
221 out << fMinY << " " << fMaxY << std::endl; <<
222 out << fMinZ << " " << fMaxZ << std::endl; <<
223 <<
224 for (unsigned int i = 0; i < fMateIDs.size() <<
225 Print(out, fMateIDs.at(i), " "); <<
226 } <<
227 <<
228 for (unsigned int i = 0; i < fValues.size(); <<
229 Print(out, fValues.at(i), " ", 6); <<
230 } <<
231 <<
232 out.close(); <<
233 } <<
234 <<
235 //....oooOO0OOooo........oooOO0OOooo........oo <<
236 void DicomPhantomZSliceHeader::ReadDataFromFil <<
237 { <<
238 std::ifstream in; <<
239 in.open(fFilename.c_str()); <<
240 <<
241 if (!in) { <<
242 G4String descript = "DicomPhantomZSliceHea <<
243 G4Exception(descript.c_str(), "", FatalExc <<
244 } <<
245 <<
246 G4int nMaterials; <<
247 in >> nMaterials; <<
248 <<
249 fMaterialNames.resize(nMaterials, ""); <<
250 for (G4int i = 0; i < nMaterials; ++i) { <<
251 G4String str1, str2; <<
252 in >> str1 >> str2; <<
253 if (!IsInteger(str1)) { <<
254 G4String descript = "String : " + str1 + <<
255 G4Exception( <<
256 "DicomPhantomZSliceHeader::ReadDataFro <<
257 formatting: missing material index", <<
258 "", FatalException, descript.c_str()); <<
259 } <<
260 G4int index = G4s2n<G4int>(str1); <<
261 if (index > nMaterials || index < 0) { <<
262 G4String descript = "Index : " + str1; <<
263 G4Exception( <<
264 "DicomPhantomZSliceHeader::ReadDataFro <<
265 bad material index", <<
266 "", FatalException, descript.c_str()); <<
267 } <<
268 fMaterialNames[index] = str2; <<
269 } <<
270 <<
271 in >> fNoVoxelsX >> fNoVoxelsY >> fNoVoxelsZ <<
272 <<
273 G4double tmpMinX, tmpMinY, tmpMinZ; <<
274 G4double tmpMaxX, tmpMaxY, tmpMaxZ; <<
275 <<
276 in >> tmpMinX >> tmpMaxX; <<
277 in >> tmpMinY >> tmpMaxY; <<
278 in >> tmpMinZ >> tmpMaxZ; <<
279 <<
280 fMinX = <<
281 (CheckConsistency(tmpMinX, fMinX, "Min X v <<
282 fMaxX = <<
283 (CheckConsistency(tmpMaxX, fMaxX, "Max X v <<
284 <<
285 fMinY = <<
286 (CheckConsistency(tmpMinY, fMinY, "Min Y v <<
287 fMaxY = <<
288 (CheckConsistency(tmpMaxY, fMaxY, "Max Y v <<
289 <<
290 fMinZ = <<
291 (CheckConsistency(tmpMinZ, fMinZ, "Min Z v <<
292 fMaxZ = <<
293 (CheckConsistency(tmpMaxZ, fMaxZ, "Max Z v <<
294 <<
295 fMateIDs.clear(); <<
296 fValues.clear(); <<
297 fMateIDs.resize(fNoVoxelsY * fNoVoxelsZ, std <<
298 fValues.resize(fNoVoxelsY * fNoVoxelsZ, std: <<
299 for (G4int k = 0; k < fNoVoxelsZ; ++k) { <<
300 for (G4int j = 0; j < fNoVoxelsY; ++j) { <<
301 for (G4int i = 0; i < fNoVoxelsX; ++i) { <<
302 G4int tmpMateID; <<
303 in >> tmpMateID; <<
304 G4int row = j * (k + 1); <<
305 fMateIDs[row][i] = tmpMateID; <<
306 } <<
307 } <<
308 } <<
309 <<
310 for (G4int k = 0; k < fNoVoxelsZ; ++k) { <<
311 for (G4int j = 0; j < fNoVoxelsY; ++j) { <<
312 for (G4int i = 0; i < fNoVoxelsX; ++i) { <<
313 G4double tmpValue; <<
314 in >> tmpValue; <<
315 G4int row = j * (k + 1); <<
316 fValues[row][i] = tmpValue; <<
317 } <<
318 } <<
319 } <<
320 <<
321 in.close(); <<
322 } <<
323 //....oooOO0OOooo........oooOO0OOooo........oo <<
324 178