| 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 // Hadrontherapy advanced example for Geant4 << 26 //
27 // See more at: https://twiki.cern.ch/twiki/bi << 27 // $Id: HadrontherapyMatrix.cc;
28 << 28 // Last modified: G.A.P.Cirrone, May 2008;
29 #include <fstream> << 29 //
30 #include <iostream> << 30 // See more at: http://geant4infn.wikispaces.com/HadrontherapyExample
31 #include <sstream> << 31 //
32 #include <iomanip> << 32 // ----------------------------------------------------------------------------
>> 33 // GEANT 4 - Hadrontherapy example
>> 34 // ----------------------------------------------------------------------------
>> 35 // Code developed by:
>> 36 //
>> 37 // G.A.P. Cirrone(a)*, F. Di Rosa(a), S. Guatelli(b), G. Russo(a)
>> 38 //
>> 39 // (a) Laboratori Nazionali del Sud
>> 40 // of the National Institute for Nuclear Physics, Catania, Italy
>> 41 // (b) National Institute for Nuclear Physics Section of Genova, genova, Italy
>> 42 //
>> 43 // * cirrone@lns.infn.it
>> 44 // ----------------------------------------------------------------------------
33 45
34 #include "HadrontherapyMatrix.hh" 46 #include "HadrontherapyMatrix.hh"
35 #include "HadrontherapyPrimaryGeneratorAction. << 47 #include "HadrontherapyAnalysisManager.hh"
36 #include "globals.hh" 48 #include "globals.hh"
37 #include "G4SystemOfUnits.hh" << 49 #include <fstream>
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 50
117 hitTrack = new G4int[numberOfVoxelAlongX*n << 51 HadrontherapyMatrix::HadrontherapyMatrix()
118 ClearHitTrack(); << 52 {
>> 53 // Number of the voxels of the phantom
>> 54 numberVoxelX = 400;
>> 55 numberVoxelY = 1;
>> 56 numberVoxelZ = 1;
>> 57
>> 58 // Create the matrix
>> 59 matrix = new G4double[numberVoxelX*numberVoxelY*numberVoxelZ];
119 } 60 }
120 61
121 ////////////////////////////////////////////// <<
122 HadrontherapyMatrix::~HadrontherapyMatrix() 62 HadrontherapyMatrix::~HadrontherapyMatrix()
123 { 63 {
124 delete[] matrix; << 64 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 } 65 }
139 <<
140 ////////////////////////////////////////////// <<
141 // Initialise the elements of the matrix to ze <<
142 <<
143 void HadrontherapyMatrix::Initialize() 66 void HadrontherapyMatrix::Initialize()
144 { << 67 {
145 // Clear ions store << 68 // Initialise the elemnts of the matrix to zero
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 69
157 ////////////////////////////////////////////// << 70 for(G4int i = 0; i < numberVoxelX; i++)
158 void HadrontherapyMatrix::PrintNuclides() <<
159 { <<
160 for (size_t i=0; i<ionStore.size(); i++) <<
161 { 71 {
162 G4cout << ionStore[i].name << G4endl; << 72 for(G4int j = 0; j < numberVoxelY; j++)
>> 73 {
>> 74 for(G4int k = 0; k < numberVoxelZ; k++)
>> 75
>> 76 matrix[(i*numberVoxelY+j)*numberVoxelZ+k] = 0.;
>> 77 }
163 } 78 }
164 } 79 }
165 80
166 ////////////////////////////////////////////// << 81 void HadrontherapyMatrix::Fill(G4int i, G4int j, G4int k,
167 // Clear Hit voxel (TrackID) markers << 82 G4double energyDeposit)
168 << 83 {
169 void HadrontherapyMatrix::ClearHitTrack() << 84 if (matrix)
170 { << 85 matrix[(i*numberVoxelY+j)*numberVoxelZ+k] += energyDeposit;
171 for(G4int i=0; i<numberOfVoxelAlongX*numbe << 86
172 } << 87 // Store the energy deposit in the matrix elemnt corresponding
173 << 88 // to the phantom voxel
174 // Return Hit status << 89 }
175 G4int* HadrontherapyMatrix::GetHitTrack(G4int << 90
176 { << 91 void HadrontherapyMatrix::TotalEnergyDeposit()
177 return &(hitTrack[Index(i,j,k)]); << 92 {
178 } << 93 // Store the information of the matrix in a ntuple and in
179 << 94 // a 1D Histogram
180 ////////////////////////////////////////////// << 95
181 // Dose methods... << 96 G4int k;
182 // Fill DOSE/fluence matrix for secondary part << 97 G4int j;
183 // If fluence parameter is true (default value << 98 G4int i;
184 // The energyDeposit parameter fill the dose m << 99
185 ////////////////////////////////////////////// << 100 if (matrix)
186 G4bool HadrontherapyMatrix::Fill(G4int trackID << 101 {
187 G4ParticleDef << 102 for(G4int l = 0; l < numberVoxelZ; l++)
188 G4int i, G4in << 103 {
189 G4double ener << 104 k = l;
190 G4bool fluenc << 105
191 { << 106 for(G4int m = 0; m < numberVoxelY; m++)
192 << 107 {
193 if ( (energyDeposit <=0. && !fluence) || ! << 108 j = m * numberVoxelZ + k;
194 << 109
195 // Get Particle Data Group particle ID << 110 for(G4int n = 0; n < numberVoxelX; n++)
196 G4int PDGencoding = particleDef -> GetPDGE << 111 {
197 PDGencoding -= PDGencoding%10; << 112 i = n* numberVoxelZ * numberVoxelY + j;
198 << 113 if(matrix[i] != 0)
199 // Search for already allocated data... << 114 {
200 for (size_t l=0; l < ionStore.size(); l++) << 115
201 { << 116 #ifdef G4ANALYSIS_USE
202 if (ionStore[l].PDGencoding == PDGenco << 117 HadrontherapyAnalysisManager* analysis =
203 { // Is it a primary or a secondary << 118 HadrontherapyAnalysisManager::getInstance();
204 << 119 analysis -> FillEnergyDeposit(n, m, k, matrix[i]);
205 if ( (trackID ==1 && ionStore[l].i << 120 analysis -> BraggPeak(n, matrix[i]);
206 { << 121 #endif
207 if (energyDeposit > 0.) << 122 }
208 << 123 }
209 ionStore[l].dose[Index(i, << 124 }
210 << 125 }
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 } 126 }
400 } 127 }
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 128