| 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 // Hadrontherapy advanced example for Geant4
27 // See more at: https://twiki.cern.ch/twiki/bi
28
29 #include <fstream>
30 #include <iostream>
31 #include <sstream>
32 #include <iomanip>
33
34 #include "HadrontherapyMatrix.hh"
35 #include "HadrontherapyPrimaryGeneratorAction.
36 #include "globals.hh"
37 #include "G4SystemOfUnits.hh"
38 #include "G4RunManager.hh"
39 #include "G4ParticleGun.hh"
40 #include "HadrontherapySteppingAction.hh"
41 #include "HadrontherapyAnalysisFileMessenger.h
42 #include "G4SystemOfUnits.hh"
43 #include <time.h>
44
45 HadrontherapyAnalysis* HadrontherapyAnalysis::
46 //////////////////////////////////////////////
47
48 HadrontherapyAnalysis::HadrontherapyAnalysis()
49 {
50 fMess = new HadrontherapyAnalysisFileMesse
51 }
52
53 //////////////////////////////////////////////
54 HadrontherapyAnalysis::~HadrontherapyAnalysis(
55 {
56 delete fMess;
57 }
58
59 //////////////////////////////////////////////
60 HadrontherapyAnalysis* HadrontherapyAnalysis::
61
62 if (instance == 0) instance = new Hadronth
63 return instance;
64 }
65
66 HadrontherapyMatrix* HadrontherapyMatrix::inst
67 G4bool HadrontherapyMatrix::secondary = false;
68
69
70 // Only return a pointer to matrix
71 HadrontherapyMatrix* HadrontherapyMatrix::GetI
72 {
73 return instance;
74 }
75
76 //////////////////////////////////////////////
77 // This STATIC method delete (!) the old matri
78 // TODO A check on the parameters is required!
79 HadrontherapyMatrix* HadrontherapyMatrix::GetI
80 {
81 if (instance) delete instance;
82 instance = new HadrontherapyMatrix(voxelX,
83 instance -> Initialize();
84 return instance;
85 }
86
87 //////////////////////////////////////////////
88 HadrontherapyMatrix::HadrontherapyMatrix(G4int
89 stdFile("Dose.out"),
90 doseUnit(gray)
91 {
92 // Number of the voxels of the phantom
93 // For Y = Z = 1 the phantom is divided in
94 // orthogonal to the beam axis
95 numberOfVoxelAlongX = voxelX;
96 numberOfVoxelAlongY = voxelY;
97 numberOfVoxelAlongZ = voxelZ;
98 massOfVoxel = mass;
99
100
101 // Create the dose matrix
102 matrix = new G4double[numberOfVoxelAlongX*
103 if (matrix)
104 {
105 G4cout << "HadrontherapyMatrix: Memory
106 numberOfVoxelAlongX*numberOf
107 " voxels has been allocated
108 }
109
110 else G4Exception("HadrontherapyMatrix::Had
111
112
113 // Hit voxel (TrackID) marker
114 // This array mark the status of voxel, if
115 // Must be initialized
116
117 hitTrack = new G4int[numberOfVoxelAlongX*n
118 ClearHitTrack();
119 }
120
121 //////////////////////////////////////////////
122 HadrontherapyMatrix::~HadrontherapyMatrix()
123 {
124 delete[] matrix;
125 delete[] hitTrack;
126 Clear();
127 }
128
129 //////////////////////////////////////////////
130 void HadrontherapyMatrix::Clear()
131 {
132 for (size_t i=0; i<ionStore.size(); i++)
133 {
134 delete[] ionStore[i].dose;
135 delete[] ionStore[i].fluence;
136 }
137 ionStore.clear();
138 }
139
140 //////////////////////////////////////////////
141 // Initialise the elements of the matrix to ze
142
143 void HadrontherapyMatrix::Initialize()
144 {
145 // Clear ions store
146 Clear();
147 // Clear dose
148 for(int i=0;i<numberOfVoxelAlongX*numberOf
149 {
150 matrix[i] = 0;
151 }
152 }
153
154 //////////////////////////////////////////////
155 // Print generated nuclides list
156
157 //////////////////////////////////////////////
158 void HadrontherapyMatrix::PrintNuclides()
159 {
160 for (size_t i=0; i<ionStore.size(); i++)
161 {
162 G4cout << ionStore[i].name << G4endl;
163 }
164 }
165
166 //////////////////////////////////////////////
167 // Clear Hit voxel (TrackID) markers
168
169 void HadrontherapyMatrix::ClearHitTrack()
170 {
171 for(G4int i=0; i<numberOfVoxelAlongX*numbe
172 }
173
174 // Return Hit status
175 G4int* HadrontherapyMatrix::GetHitTrack(G4int
176 {
177 return &(hitTrack[Index(i,j,k)]);
178 }
179
180 //////////////////////////////////////////////
181 // Dose methods...
182 // Fill DOSE/fluence matrix for secondary part
183 // If fluence parameter is true (default value
184 // The energyDeposit parameter fill the dose m
185 //////////////////////////////////////////////
186 G4bool HadrontherapyMatrix::Fill(G4int trackID
187 G4ParticleDef
188 G4int i, G4in
189 G4double ener
190 G4bool fluenc
191 {
192
193 if ( (energyDeposit <=0. && !fluence) || !
194
195 // Get Particle Data Group particle ID
196 G4int PDGencoding = particleDef -> GetPDGE
197 PDGencoding -= PDGencoding%10;
198
199 // Search for already allocated data...
200 for (size_t l=0; l < ionStore.size(); l++)
201 {
202 if (ionStore[l].PDGencoding == PDGenco
203 { // Is it a primary or a secondary
204
205 if ( (trackID ==1 && ionStore[l].i
206 {
207 if (energyDeposit > 0.)
208
209 ionStore[l].dose[Index(i,
210
211 // Fill a matrix per each ion
212
213 if (fluence) ionStore[l].fluen
214 return true;
215 }
216 }
217 }
218 G4int Z = particleDef-> GetAtomicNumber();
219 G4int A = particleDef-> GetAtomicMass();
220 G4String fullName = particleDef -> GetPart
221 G4String name = fullName.substr (0, fullNa
222
223 // Let's put a new particle in our store..
224 ion newIon =
225 {
226 (trackID == 1) ? true:false,
227 PDGencoding,
228 name,
229 name.length(),
230 Z,
231 A,
232 new G4double[numberOfVoxelAlongX * num
233 new unsigned int[numberOfVoxelAlongX *
234 };
235
236
237 // Initialize data
238 if (newIon.dose && newIon.fluence)
239 {
240 for(G4int q=0; q<numberOfVoxelAlongX*n
241 {
242 newIon.dose[q] = 0.;
243 newIon.fluence[q] = 0;
244 }
245
246 if (energyDeposit > 0.) newIon.dose[In
247 if (fluence) newIon.fluence[Index(i, j
248
249 ionStore.push_back(newIon);
250 return true;
251 }
252
253 else // XXX Out of memory! XXX
254
255 {
256 return false;
257 }
258 }
259
260 //////////////////////////////////////////////
261 //////////////////////////////////////////////
262 // Methods to store data to filenames...
263 //////////////////////////////////////////////
264 //////////////////////////////////////////////
265 //
266 // General method to store matrix data to file
267 void HadrontherapyMatrix::StoreMatrix(G4String
268 {
269 if (data)
270 {
271 ofs.open(file, std::ios::out);
272 if (ofs.is_open())
273 {
274 for(G4int i = 0; i < numberOfVoxel
275 for(G4int j = 0; j < numberOfV
276 for(G4int k = 0; k < numbe
277 {
278 G4int n = Index(i, j,
279
280 if (psize == sizeof(un
281 {
282 unsigned int* pdat
283
284 if (pdata[n])
285
286 ofs << i << '\
287 }
288
289 else if (psize == size
290
291 {
292 G4double* pdata =
293 if (pdata[n]) ofs
294 }
295 }
296 ofs.close();
297 }
298 }
299 }
300
301 //////////////////////////////////////////////
302 // Store fluence per single ion in multiple fi
303 void HadrontherapyMatrix::StoreFluenceData()
304 {
305 for (size_t i=0; i < ionStore.size(); i++)
306 StoreMatrix(ionStore[i].name + "_Fluen
307 }
308 }
309
310 //////////////////////////////////////////////
311 // Store dose per single ion in multiple files
312 void HadrontherapyMatrix::StoreDoseData()
313 {
314
315 for (size_t i=0; i < ionStore.size(); i++)
316 StoreMatrix(ionStore[i].name + "_Dose.
317 }
318 }
319
320 //////////////////////////////////////////////
321 // Store dose into a single file
322 // or in histograms. Please, note that this fu
323 // messenger commands
324 // defined in the HadrontherapyAnalysisFileMes
325 void HadrontherapyMatrix::StoreDoseFluenceAsci
326 {
327 #define width 15L
328 filename = (file=="") ? stdFile:file;
329
330 // Sort like periodic table
331
332 std::sort(ionStore.begin(), ionStore.end()
333 G4cout << "Dose is being written to " << f
334 ofs.open(filename, std::ios::out);
335
336 if (ofs.is_open())
337 {
338 // Write the voxels index and the list
339 //ofs << std::setprecision(6) << std::
340 ofs << "i" << '\t' << "j" << '\t' << "
341 G4cout << "i" << '\t' << "j" << '\t' <
342
343 // Total dose
344 ofs <<'\t' <<"Dose(Gy)";
345 //ofs << std::setw(width) << "Dose(Gy)
346 G4cout << '\t' << "Dose(Gy)";
347
348 G4String fluence = "_f";
349 if (secondary)
350 {
351 for (size_t l=0; l < ionStore.size
352 {
353 G4String a = (ionStore[l].isPr
354
355 // ofs << std::setw(width) << i
356 // std::setw(width) << i
357
358 ofs << '\t' << ionStore[l].nam
359 '\t' << ionStore[l].
360
361 G4cout << '\t' << ionStore[l].
362 '\t' << ionStore[l].
363
364 }
365 //ofs << G4endl;
366 }
367
368 // Write data
369 for(G4int i = 0; i < numberOfVoxelAlon
370 for(G4int j = 0; j < numberOfVoxel
371 for(G4int k = 0; k < numberOfV
372 {
373 G4int n = Index(i, j, k);
374
375 if (matrix[n])
376 {
377 ofs << G4endl;
378 ofs << i << '\t' << j
379
380 // Total dose
381 //ofs << std::setw(wid
382 ofs << (matrix[n]/mass
383
384 if (secondary)
385 {
386 for (size_t l=0; l
387 {
388 // Fill ASCII
389 //ofs << std::
390 // std::
391
392 ofs << '\t' <
393 '\t' << i
394 }
395 }
396 }
397 }
398 ofs.close();
399 }
400 }
401 //////////////////////////////////////////////
402 void HadrontherapyMatrix::Fill(G4int i, G4int
403 G4double energy
404 {
405 if (matrix)
406 matrix[Index(i,j,k)] += energyDeposit;
407
408 // Store the energy deposit in the matrix
409 // to the phantom voxel
410 }
411
412
413
414
415