1 | // |
---|
2 | // ******************************************************************** |
---|
3 | // * License and Disclaimer * |
---|
4 | // * * |
---|
5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
7 | // * conditions of the Geant4 Software License, included in the file * |
---|
8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
9 | // * include a list of copyright holders. * |
---|
10 | // * * |
---|
11 | // * Neither the authors of this software system, nor their employing * |
---|
12 | // * institutes,nor the agencies providing financial support for this * |
---|
13 | // * work make any representation or warranty, express or implied, * |
---|
14 | // * regarding this software system or assume any liability for its * |
---|
15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
16 | // * for the full disclaimer and the limitation of liability. * |
---|
17 | // * * |
---|
18 | // * This code implementation is the result of the scientific and * |
---|
19 | // * technical work of the GEANT4 collaboration. * |
---|
20 | // * By using, copying, modifying or distributing the software (or * |
---|
21 | // * any work based on the software) you agree to acknowledge its * |
---|
22 | // * use in resulting scientific publications, and indicate your * |
---|
23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
24 | // ******************************************************************** |
---|
25 | // |
---|
26 | // $Id: HadrontherapyMatrix.cc; |
---|
27 | // See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy// |
---|
28 | |
---|
29 | #include "HadrontherapyMatrix.hh" |
---|
30 | #include "HadrontherapyAnalysisManager.hh" |
---|
31 | #include "G4RunManager.hh" |
---|
32 | #include "HadrontherapyPrimaryGeneratorAction.hh" |
---|
33 | #include "G4ParticleGun.hh" |
---|
34 | |
---|
35 | #include <fstream> |
---|
36 | #include <unistd.h> |
---|
37 | #include <iostream> |
---|
38 | #include <sstream> |
---|
39 | #include <iomanip> |
---|
40 | #include "globals.hh" |
---|
41 | |
---|
42 | // Units definition: CLHEP/Units/SystemOfUnits.h |
---|
43 | // |
---|
44 | HadrontherapyMatrix* HadrontherapyMatrix::instance = NULL; |
---|
45 | G4bool HadrontherapyMatrix::secondaries = false; |
---|
46 | // Only return a pointer to matrix |
---|
47 | HadrontherapyMatrix* HadrontherapyMatrix::GetInstance() |
---|
48 | { |
---|
49 | return instance; |
---|
50 | } |
---|
51 | // This STATIC method delete (!) the old matrix and rewrite a new object returning a pointer to it |
---|
52 | // TODO A check on the parameters is required! |
---|
53 | HadrontherapyMatrix* HadrontherapyMatrix::GetInstance(G4int voxelX, G4int voxelY, G4int voxelZ, G4double mass) |
---|
54 | { |
---|
55 | if (instance) delete instance; |
---|
56 | instance = new HadrontherapyMatrix(voxelX, voxelY, voxelZ, mass); |
---|
57 | instance -> Initialize();// XXX |
---|
58 | return instance; |
---|
59 | } |
---|
60 | HadrontherapyMatrix::HadrontherapyMatrix(G4int voxelX, G4int voxelY, G4int voxelZ, G4double mass): |
---|
61 | filename("Dose.out"), |
---|
62 | doseUnit(MeV/g) |
---|
63 | { |
---|
64 | // Number of the voxels of the phantom |
---|
65 | // For Y = Z = 1 the phantom is divided in slices (and not in voxels) |
---|
66 | // orthogonal to the beam axis |
---|
67 | numberOfVoxelAlongX = voxelX; |
---|
68 | numberOfVoxelAlongY = voxelY; |
---|
69 | numberOfVoxelAlongZ = voxelZ; |
---|
70 | massOfVoxel = mass; |
---|
71 | // Create the dose matrix |
---|
72 | matrix = new G4double[numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ]; |
---|
73 | if (matrix) |
---|
74 | { |
---|
75 | G4cout << "HadrontherapyMatrix: Memory space to store physical dose into " << |
---|
76 | numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ << |
---|
77 | " voxels has been allocated " << G4endl; |
---|
78 | } |
---|
79 | else G4Exception(" HadrontherapyMatrix::HadrontherapyMatrix. Can't allocate memory to store physical dose!"); |
---|
80 | // Hit voxel (TrackID) marker |
---|
81 | // This array mark the status of voxel, if a hit occur, with the trackID of the particle |
---|
82 | // Must be initialized |
---|
83 | hitTrack = new G4int[numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ]; |
---|
84 | ClearHitTrack(); |
---|
85 | } |
---|
86 | |
---|
87 | ///////////////////////////////////////////////////////////////////////////// |
---|
88 | HadrontherapyMatrix::~HadrontherapyMatrix() |
---|
89 | { |
---|
90 | delete[] matrix; |
---|
91 | delete[] hitTrack; |
---|
92 | // clear fluences/dose data |
---|
93 | Clear(); |
---|
94 | } |
---|
95 | |
---|
96 | ///////////////////////////////////////////////////////////////////////////// |
---|
97 | void HadrontherapyMatrix::Clear() |
---|
98 | { |
---|
99 | for (size_t i=0; i<ionStore.size(); i++) |
---|
100 | { |
---|
101 | delete[] ionStore[i].dose; |
---|
102 | delete[] ionStore[i].fluence; |
---|
103 | } |
---|
104 | ionStore.clear(); |
---|
105 | } |
---|
106 | |
---|
107 | ///////////////////////////////////////////////////////////////////////////// |
---|
108 | // Initialise the elements of the matrix to zero |
---|
109 | void HadrontherapyMatrix::Initialize() |
---|
110 | { |
---|
111 | Clear(); |
---|
112 | for(int i=0;i<numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ;i++) |
---|
113 | { |
---|
114 | matrix[i] = 0; |
---|
115 | } |
---|
116 | } |
---|
117 | ///////////////////////////////////////////////////////////////////////////// |
---|
118 | ///////////////////////////////////////////////////////////////////////////// |
---|
119 | // Print generated nuclides list |
---|
120 | void HadrontherapyMatrix::PrintNuclides() |
---|
121 | { |
---|
122 | for (size_t i=0; i<ionStore.size(); i++) |
---|
123 | { |
---|
124 | G4cout << ionStore[i].name << G4endl; |
---|
125 | } |
---|
126 | } |
---|
127 | ///////////////////////////////////////////////////////////////////////////// |
---|
128 | // Clear Hit voxel (TrackID) markers |
---|
129 | void HadrontherapyMatrix::ClearHitTrack() |
---|
130 | { |
---|
131 | for(G4int i=0; i<numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ; i++) hitTrack[i] = 0; |
---|
132 | } |
---|
133 | // Return Hit status |
---|
134 | G4int* HadrontherapyMatrix::GetHitTrack(G4int i, G4int j, G4int k) |
---|
135 | { |
---|
136 | return &(hitTrack[Index(i,j,k)]); |
---|
137 | } |
---|
138 | ///////////////////////////////////////////////////////////////////////////// |
---|
139 | // Dose methods... |
---|
140 | // Fill DOSE/fluence matrix for particle/ion: |
---|
141 | // If fluence parameter is true (default value is FALSE) then fluence at voxel (i, j, k) is increased. |
---|
142 | // The energyDeposit parameter fill the dose matrix for voxel (i,j,k) |
---|
143 | ///////////////////////////////////////////////////////////////////////////// |
---|
144 | |
---|
145 | G4bool HadrontherapyMatrix::Fill(G4int trackID, |
---|
146 | G4ParticleDefinition* particleDef, |
---|
147 | G4int i, G4int j, G4int k, |
---|
148 | G4double energyDeposit, |
---|
149 | G4bool fluence) |
---|
150 | { |
---|
151 | if (energyDeposit <=0. && !fluence || !secondaries) return false; |
---|
152 | // Get Particle Data Group particle ID |
---|
153 | G4int PDGencoding = particleDef -> GetPDGEncoding(); |
---|
154 | PDGencoding -= PDGencoding%10; |
---|
155 | |
---|
156 | // Search for already allocated data... |
---|
157 | for (size_t l=0; l < ionStore.size(); l++) |
---|
158 | { |
---|
159 | if (ionStore[l].PDGencoding == PDGencoding ) |
---|
160 | { // Is it a primary or a secondary particle? |
---|
161 | if ( trackID ==1 && ionStore[l].isPrimary || trackID !=1 && !ionStore[l].isPrimary) |
---|
162 | { |
---|
163 | if (energyDeposit > 0.) ionStore[l].dose[Index(i, j, k)] += energyDeposit/massOfVoxel; |
---|
164 | |
---|
165 | // Fill a matrix per each ion with the fluence |
---|
166 | if (fluence) ionStore[l].fluence[Index(i, j, k)]++; |
---|
167 | return true; |
---|
168 | } |
---|
169 | } |
---|
170 | } |
---|
171 | |
---|
172 | G4int Z = particleDef-> GetAtomicNumber(); |
---|
173 | G4int A = particleDef-> GetAtomicMass(); |
---|
174 | |
---|
175 | G4String fullName = particleDef -> GetParticleName(); |
---|
176 | G4String name = fullName.substr (0, fullName.find("[") ); // cut excitation energy |
---|
177 | // Let's put a new particle in our store... |
---|
178 | ion newIon = |
---|
179 | { |
---|
180 | (trackID == 1) ? true:false, |
---|
181 | PDGencoding, |
---|
182 | name, |
---|
183 | name.length(), |
---|
184 | Z, |
---|
185 | A, |
---|
186 | new G4double[numberOfVoxelAlongX * numberOfVoxelAlongY * numberOfVoxelAlongZ], |
---|
187 | new unsigned int[numberOfVoxelAlongX * numberOfVoxelAlongY * numberOfVoxelAlongZ] |
---|
188 | }; |
---|
189 | // Initialize data |
---|
190 | if (newIon.dose && newIon.fluence) |
---|
191 | { |
---|
192 | for(G4int m=0; m<numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ; m++) |
---|
193 | { |
---|
194 | newIon.dose[m] = 0.; |
---|
195 | newIon.fluence[m] = 0; |
---|
196 | } |
---|
197 | if (energyDeposit > 0.) newIon.dose[Index(i, j, k)] += energyDeposit/massOfVoxel; |
---|
198 | if (fluence) newIon.fluence[Index(i, j, k)]++; |
---|
199 | |
---|
200 | ionStore.push_back(newIon); |
---|
201 | |
---|
202 | // TODO Put some verbosity check |
---|
203 | /* |
---|
204 | G4cout << "Memory space to store the DOSE/FLUENCE into " << |
---|
205 | numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ << |
---|
206 | " voxels has been allocated for the nuclide " << newIon.name << |
---|
207 | " (Z = " << Z << ", A = " << A << ")" << G4endl ; |
---|
208 | */ |
---|
209 | return true; |
---|
210 | } |
---|
211 | else // XXX Out of memory! XXX |
---|
212 | { |
---|
213 | return false; |
---|
214 | } |
---|
215 | |
---|
216 | } |
---|
217 | |
---|
218 | ///////////////////////////////////////////////////////////////////////////// |
---|
219 | ///////////////////////////////////////////////////////////////////////////// |
---|
220 | // Methods to store data to filenames... |
---|
221 | //////////////////////////////////////////////////////////////////////////// |
---|
222 | //////////////////////////////////////////////////////////////////////////// |
---|
223 | // |
---|
224 | // General method to store matrix data to filename |
---|
225 | void HadrontherapyMatrix::StoreMatrix(G4String file, void* data, size_t psize) |
---|
226 | { |
---|
227 | if (data) |
---|
228 | { |
---|
229 | ofs.open(file, std::ios::out); |
---|
230 | if (ofs.is_open()) |
---|
231 | { |
---|
232 | for(G4int i = 0; i < numberOfVoxelAlongX; i++) |
---|
233 | for(G4int j = 0; j < numberOfVoxelAlongY; j++) |
---|
234 | for(G4int k = 0; k < numberOfVoxelAlongZ; k++) |
---|
235 | { |
---|
236 | G4int n = Index(i, j, k); |
---|
237 | // Check for data type: u_int, G4double, XXX |
---|
238 | if (psize == sizeof(unsigned int)) |
---|
239 | { |
---|
240 | unsigned int* pdata = (unsigned int*)data; |
---|
241 | if (pdata[n]) ofs << i << '\t' << j << '\t' << |
---|
242 | k << '\t' << pdata[n] << G4endl; |
---|
243 | } |
---|
244 | else if (psize == sizeof(G4double)) |
---|
245 | { |
---|
246 | G4double* pdata = (G4double*)data; |
---|
247 | if (pdata[n]) ofs << i << '\t' << j << '\t' << |
---|
248 | k << '\t' << pdata[n] << G4endl; |
---|
249 | } |
---|
250 | } |
---|
251 | ofs.close(); |
---|
252 | } |
---|
253 | } |
---|
254 | } |
---|
255 | |
---|
256 | // Store fluence per single ion in multiple files |
---|
257 | void HadrontherapyMatrix::StoreFluenceData() |
---|
258 | { |
---|
259 | for (size_t i=0; i < ionStore.size(); i++){ |
---|
260 | StoreMatrix(ionStore[i].name + "_Fluence.out", ionStore[i].fluence, sizeof(unsigned int)); |
---|
261 | } |
---|
262 | } |
---|
263 | // Store dose per single ion in multiple files |
---|
264 | void HadrontherapyMatrix::StoreDoseData() |
---|
265 | { |
---|
266 | |
---|
267 | for (size_t i=0; i < ionStore.size(); i++){ |
---|
268 | StoreMatrix(ionStore[i].name + "_Dose.out", ionStore[i].dose, sizeof(G4double)); |
---|
269 | } |
---|
270 | } |
---|
271 | ///////////////////////////////////////////////////////////////////////// |
---|
272 | // Store dose for all ions into a single file and into ntuples. |
---|
273 | // Please note that this function is called via messenger commands |
---|
274 | // defined in the HadrontherapyAnalysisFileMessenger.cc class file |
---|
275 | void HadrontherapyMatrix::StoreDoseFluenceAscii() |
---|
276 | { |
---|
277 | // Sort like periodic table |
---|
278 | std::sort(ionStore.begin(), ionStore.end()); |
---|
279 | #define width 15L |
---|
280 | ofs.open(filename, std::ios::out); |
---|
281 | if (ofs.is_open()) |
---|
282 | { |
---|
283 | |
---|
284 | // |
---|
285 | // Write the voxels index and the list of particles/ions |
---|
286 | ofs << std::setprecision(6) << std::left << |
---|
287 | "i\tj\tk\t"; |
---|
288 | /* |
---|
289 | G4RunManager *runManager = G4RunManager::GetRunManager(); |
---|
290 | HadrontherapyPrimaryGeneratorAction *pPGA = (HadrontherapyPrimaryGeneratorAction*)runManager -> GetUserPrimaryGeneratorAction(); |
---|
291 | G4String name = pPGA -> GetParticleGun() -> GetParticleDefinition() -> GetParticleName(); |
---|
292 | */ |
---|
293 | // Total dose |
---|
294 | ofs << std::setw(width) << "Dose"; |
---|
295 | |
---|
296 | for (size_t l=0; l < ionStore.size(); l++) |
---|
297 | { |
---|
298 | G4String a = (ionStore[l].isPrimary) ? "_1":""; // is it a primary? |
---|
299 | ofs << std::setw(width) << ionStore[l].name + a << |
---|
300 | std::setw(width) << ionStore[l].name + a; |
---|
301 | } |
---|
302 | ofs << G4endl; |
---|
303 | |
---|
304 | /* |
---|
305 | // PDGencoding |
---|
306 | ofs << std::setprecision(6) << std::left << |
---|
307 | "0\t0\t0\t"; |
---|
308 | |
---|
309 | // Total dose |
---|
310 | ofs << std::setw(width) << '0'; |
---|
311 | for (size_t l=0; l < ionStore.size(); l++) |
---|
312 | { |
---|
313 | ofs << std::setw(width) << ionStore[l].PDGencoding << |
---|
314 | std::setw(width) << ionStore[l].PDGencoding; |
---|
315 | } |
---|
316 | ofs << G4endl; |
---|
317 | */ |
---|
318 | // Write data |
---|
319 | for(G4int i = 0; i < numberOfVoxelAlongX; i++) |
---|
320 | for(G4int j = 0; j < numberOfVoxelAlongY; j++) |
---|
321 | for(G4int k = 0; k < numberOfVoxelAlongZ; k++) |
---|
322 | { |
---|
323 | G4int n = Index(i, j, k); |
---|
324 | // Write only not identically null data lines |
---|
325 | if (matrix[n]) |
---|
326 | { |
---|
327 | ofs << G4endl; |
---|
328 | ofs << i << '\t' << j << '\t' << k << '\t'; |
---|
329 | // Total dose |
---|
330 | ofs << std::setw(width) << matrix[n]/(doseUnit); |
---|
331 | { |
---|
332 | for (size_t l=0; l < ionStore.size(); l++) |
---|
333 | { |
---|
334 | // Fill ASCII file rows |
---|
335 | ofs << std::setw(width) << ionStore[l].dose[n]/(doseUnit) << |
---|
336 | std::setw(width) << ionStore[l].fluence[n]; |
---|
337 | } |
---|
338 | } |
---|
339 | } |
---|
340 | } |
---|
341 | ofs.close(); |
---|
342 | } |
---|
343 | } |
---|
344 | ///////////////////////////////////////////////////////////////////////////// |
---|
345 | |
---|
346 | #ifdef G4ANALYSIS_USE_ROOT |
---|
347 | void HadrontherapyMatrix::StoreDoseFluenceRoot() |
---|
348 | { |
---|
349 | HadrontherapyAnalysisManager* analysis = HadrontherapyAnalysisManager::GetInstance(); |
---|
350 | for(G4int i = 0; i < numberOfVoxelAlongX; i++) |
---|
351 | for(G4int j = 0; j < numberOfVoxelAlongY; j++) |
---|
352 | for(G4int k = 0; k < numberOfVoxelAlongZ; k++) |
---|
353 | { |
---|
354 | G4int n = Index(i, j, k); |
---|
355 | for (size_t l=0; l < ionStore.size(); l++) |
---|
356 | |
---|
357 | { |
---|
358 | // Do the same work for .root file: fill dose/fluence ntuple |
---|
359 | analysis -> FillVoxelFragmentTuple( i, j, k, |
---|
360 | ionStore[l].A, |
---|
361 | ionStore[l].Z, |
---|
362 | ionStore[l].dose[n]/(doseUnit), |
---|
363 | ionStore[l].fluence[n] ); |
---|
364 | |
---|
365 | |
---|
366 | } |
---|
367 | } |
---|
368 | } |
---|
369 | #endif |
---|
370 | |
---|
371 | void HadrontherapyMatrix::Fill(G4int i, G4int j, G4int k, |
---|
372 | G4double energyDeposit) |
---|
373 | { |
---|
374 | if (matrix) |
---|
375 | matrix[Index(i,j,k)] += energyDeposit; |
---|
376 | |
---|
377 | // Store the energy deposit in the matrix element corresponding |
---|
378 | // to the phantom voxel |
---|
379 | } |
---|
380 | void HadrontherapyMatrix::TotalEnergyDeposit() |
---|
381 | { |
---|
382 | // Convert energy deposited to dose. |
---|
383 | // Store the information of the matrix in a ntuple and in |
---|
384 | // a 1D Histogram |
---|
385 | |
---|
386 | #ifdef G4ANALYSIS_USE_ROOT |
---|
387 | HadrontherapyAnalysisManager* analysis = HadrontherapyAnalysisManager::GetInstance(); |
---|
388 | #endif |
---|
389 | if (matrix) |
---|
390 | { |
---|
391 | for(G4int i = 0; i < numberOfVoxelAlongX; i++) |
---|
392 | for(G4int j = 0; j < numberOfVoxelAlongY; j++) |
---|
393 | for(G4int k = 0; k < numberOfVoxelAlongZ; k++) |
---|
394 | { |
---|
395 | G4int n = Index(i,j,k); |
---|
396 | matrix[n] = matrix[n] / massOfVoxel; |
---|
397 | #ifdef G4ANALYSIS_USE_ROOT |
---|
398 | analysis -> FillEnergyDeposit(i, j, k, matrix[n]/doseUnit); |
---|
399 | analysis -> BraggPeak(i, matrix[n]/doseUnit); |
---|
400 | #endif |
---|
401 | } |
---|
402 | } |
---|
403 | } |
---|
404 | |
---|