| 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 // $Id: HadrontherapyPhantomSD.cc; May 2005
27 // See more at: https://twiki.cern.ch/twiki/bi << 27 // ----------------------------------------------------------------------------
28 << 28 // GEANT 4 - Hadrontherapy example
29 #include <fstream> << 29 // ----------------------------------------------------------------------------
30 #include <iostream> << 30 // Code developed by:
31 #include <sstream> << 31 //
32 #include <iomanip> << 32 // G.A.P. Cirrone(a)*, F. Di Rosa(a), S. Guatelli(b), G. Russo(a)
>> 33 //
>> 34 // (a) Laboratori Nazionali del Sud
>> 35 // of the National Institute for Nuclear Physics, Catania, Italy
>> 36 // (b) National Institute for Nuclear Physics Section of Genova, genova, Italy
>> 37 //
>> 38 // * cirrone@lns.infn.it
>> 39 // ----------------------------------------------------------------------------
33 40
34 #include "HadrontherapyMatrix.hh" 41 #include "HadrontherapyMatrix.hh"
35 #include "HadrontherapyPrimaryGeneratorAction. << 42 #include "HadrontherapyAnalysisManager.hh"
36 #include "globals.hh" 43 #include "globals.hh"
37 #include "G4SystemOfUnits.hh" << 44 #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 45
117 hitTrack = new G4int[numberOfVoxelAlongX*n << 46 HadrontherapyMatrix::HadrontherapyMatrix()
118 ClearHitTrack(); << 47 {
>> 48 // Number of the voxels of the phantom
>> 49 numberVoxelX = 200;
>> 50 numberVoxelY = 200;
>> 51 numberVoxelZ = 200;
>> 52
>> 53 // Create the matrix
>> 54 matrix = new G4double[numberVoxelX*numberVoxelY*numberVoxelZ];
119 } 55 }
120 56
121 ////////////////////////////////////////////// <<
122 HadrontherapyMatrix::~HadrontherapyMatrix() 57 HadrontherapyMatrix::~HadrontherapyMatrix()
123 { 58 {
124 delete[] matrix; << 59 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 } 60 }
139 <<
140 ////////////////////////////////////////////// <<
141 // Initialise the elements of the matrix to ze <<
142 <<
143 void HadrontherapyMatrix::Initialize() 61 void HadrontherapyMatrix::Initialize()
144 { << 62 {
145 // Clear ions store << 63 // 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 <<
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 64
199 // Search for already allocated data... << 65 for(G4int i = 0; i < numberVoxelX; i++)
200 for (size_t l=0; l < ionStore.size(); l++) <<
201 { 66 {
202 if (ionStore[l].PDGencoding == PDGenco << 67 for(G4int j = 0; j < numberVoxelY; j++)
203 { // Is it a primary or a secondary << 68 {
204 << 69 for(G4int k = 0; k < numberVoxelZ; k++)
205 if ( (trackID ==1 && ionStore[l].i << 70
206 { << 71 matrix[(i*numberVoxelY+j)*numberVoxelZ+k] = 0.;
207 if (energyDeposit > 0.) << 72 }
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 } 73 }
258 } 74 }
259 75
260 ////////////////////////////////////////////// << 76 void HadrontherapyMatrix::Fill(G4int i, G4int j, G4int k,
261 ////////////////////////////////////////////// << 77 G4double energyDeposit)
262 // Methods to store data to filenames... << 78 {
263 ////////////////////////////////////////////// << 79 if (matrix)
264 ////////////////////////////////////////////// << 80 matrix[(i*numberVoxelY+j)*numberVoxelZ+k] += energyDeposit;
265 // << 81
266 // General method to store matrix data to file << 82 // Store the energy deposit in the matrix elemnt corresponding
267 void HadrontherapyMatrix::StoreMatrix(G4String << 83 // to the phantom voxel
268 { << 84 }
269 if (data) << 85
270 { << 86 void HadrontherapyMatrix::TotalEnergyDeposit()
271 ofs.open(file, std::ios::out); << 87 {
272 if (ofs.is_open()) << 88 // Store the information of the matrix in a ntuple and in
273 { << 89 // a 1D Histogram
274 for(G4int i = 0; i < numberOfVoxel << 90
275 for(G4int j = 0; j < numberOfV << 91 G4int k;
276 for(G4int k = 0; k < numbe << 92 G4int j;
277 { << 93 G4int i;
278 G4int n = Index(i, j, << 94
279 << 95 if (matrix)
280 if (psize == sizeof(un << 96 {
281 { << 97 for(G4int l = 0; l < numberVoxelZ; l++)
282 unsigned int* pdat << 98 {
283 << 99 k = l;
284 if (pdata[n]) << 100
285 << 101 for(G4int m = 0; m < numberVoxelY; m++)
286 ofs << i << '\ << 102 {
287 } << 103 j = m * numberVoxelZ + k;
288 << 104
289 else if (psize == size << 105 for(G4int n = 0; n < numberVoxelX; n++)
290 << 106 {
291 { << 107 i = n* numberVoxelZ * numberVoxelY + j;
292 G4double* pdata = << 108 if(matrix[i] != 0)
293 if (pdata[n]) ofs << 109 {
294 } << 110
295 } << 111 #ifdef G4ANALYSIS_USE
296 ofs.close(); << 112 HadrontherapyAnalysisManager* analysis =
297 } << 113 HadrontherapyAnalysisManager::getInstance();
>> 114 analysis -> FillEnergyDeposit(n, m, k, matrix[i]);
>> 115 analysis -> BraggPeak(n, matrix[i]);
>> 116 #endif
>> 117 }
>> 118 }
>> 119 }
>> 120 }
298 } 121 }
299 } 122 }
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 123