ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
HadrontherapyMatrix.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file HadrontherapyMatrix.cc
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 // Hadrontherapy advanced example for Geant4
27 // See more at: https://twiki.cern.ch/twiki/bin/view/Geant4/AdvancedExamplesHadrontherapy
28 
29 #include <fstream>
30 #include <iostream>
31 #include <sstream>
32 #include <iomanip>
33 
34 #include "HadrontherapyMatrix.hh"
36 #include "globals.hh"
37 #include "G4SystemOfUnits.hh"
38 #include "G4RunManager.hh"
39 #include "G4ParticleGun.hh"
42 #include "G4SystemOfUnits.hh"
43 #include <time.h>
44 
47 
49 {
51 }
52 
55 {
56  delete fMess;
57 }
58 
61 
62  if (instance == 0) instance = new HadrontherapyAnalysis;
63  return instance;
64 }
65 
68 
69 
70 // Only return a pointer to matrix
72 {
73  return instance;
74 }
75 
77 // This STATIC method delete (!) the old matrix and rewrite a new object returning a pointer to it
78 // TODO A check on the parameters is required!
80 {
81  if (instance) delete instance;
82  instance = new HadrontherapyMatrix(voxelX, voxelY, voxelZ, mass);
83  instance -> Initialize();
84  return instance;
85 }
86 
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 slices (and not in voxels)
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
103  if (matrix)
104  {
105  G4cout << "HadrontherapyMatrix: Memory space to store physical dose into " <<
107  " voxels has been allocated " << G4endl;
108  }
109 
110 
111  else G4Exception("HadrontherapyMatrix::HadrontherapyMatrix()", "Hadrontherapy0005", FatalException, "Can't allocate memory to store physical dose!");
112 
113 
114  // Hit voxel (TrackID) marker
115  // This array mark the status of voxel, if a hit occur, with the trackID of the particle
116  // Must be initialized
117 
119  ClearHitTrack();
120 }
121 
124 {
125  delete[] matrix;
126  delete[] hitTrack;
127  Clear();
128 }
129 
132 {
133  for (size_t i=0; i<ionStore.size(); i++)
134  {
135  delete[] ionStore[i].dose;
136  delete[] ionStore[i].fluence;
137  }
138  ionStore.clear();
139 }
140 
142 // Initialise the elements of the matrix to zero
143 
145 {
146  // Clear ions store
147  Clear();
148  // Clear dose
150  {
151  matrix[i] = 0;
152  }
153 }
154 
156 // Print generated nuclides list
157 
160 {
161  for (size_t i=0; i<ionStore.size(); i++)
162  {
163  G4cout << ionStore[i].name << G4endl;
164  }
165 }
166 
168 // Clear Hit voxel (TrackID) markers
169 
171 {
173 }
174 
175 
176 // Return Hit status
178 {
179  return &(hitTrack[Index(i,j,k)]);
180 }
181 
182 
183 
185 // Dose methods...
186 // Fill DOSE/fluence matrix for secondary particles:
187 // If fluence parameter is true (default value is FALSE) then fluence at voxel (i, j, k) is increased.
188 // The energyDeposit parameter fill the dose matrix for voxel (i,j,k)
191  G4ParticleDefinition* particleDef,
192  G4int i, G4int j, G4int k,
193  G4double energyDeposit,
194  G4bool fluence)
195 {
196 
197  if ( (energyDeposit <=0. && !fluence) || !secondary) return false;
198 
199  // Get Particle Data Group particle ID
200  G4int PDGencoding = particleDef -> GetPDGEncoding();
201  PDGencoding -= PDGencoding%10;
202 
203  // Search for already allocated data...
204  for (size_t l=0; l < ionStore.size(); l++)
205  {
206  if (ionStore[l].PDGencoding == PDGencoding )
207  { // Is it a primary or a secondary particle?
208 
209  if ( (trackID ==1 && ionStore[l].isPrimary) || (trackID !=1 && !ionStore[l].isPrimary))
210  {
211  if (energyDeposit > 0.)
212 
213  ionStore[l].dose[Index(i, j, k)] += energyDeposit;
214 
215  // Fill a matrix per each ion with the fluence
216 
217  if (fluence) ionStore[l].fluence[Index(i, j, k)]++;
218  return true;
219  }
220  }
221  }
222 
223  G4int Z = particleDef-> GetAtomicNumber();
224  G4int A = particleDef-> GetAtomicMass();
225  G4String fullName = particleDef -> GetParticleName();
226  G4String name = fullName.substr (0, fullName.find("[") ); // cut excitation energy
227 
228  // Let's put a new particle in our store...
229 
230 
231  ion newIon =
232  {
233  (trackID == 1) ? true:false,
234  PDGencoding,
235  name,
236  name.length(),
237  Z,
238  A,
241  };
242 
243 
244  // Initialize data
245  if (newIon.dose && newIon.fluence)
246  {
248  {
249  newIon.dose[q] = 0.;
250  newIon.fluence[q] = 0;
251  }
252 
253  if (energyDeposit > 0.) newIon.dose[Index(i, j, k)] += energyDeposit;
254  if (fluence) newIon.fluence[Index(i, j, k)]++;
255 
256  ionStore.push_back(newIon);
257  return true;
258  }
259 
260  else // XXX Out of memory! XXX
261 
262  {
263  return false;
264  }
265 
266 }
267 
270 // Methods to store data to filenames...
273 //
274 // General method to store matrix data to filename
276 {
277  if (data)
278  {
279  ofs.open(file, std::ios::out);
280  if (ofs.is_open())
281  {
282  for(G4int i = 0; i < numberOfVoxelAlongX; i++)
283  for(G4int j = 0; j < numberOfVoxelAlongY; j++)
284  for(G4int k = 0; k < numberOfVoxelAlongZ; k++)
285  {
286  G4int n = Index(i, j, k);
287 
288  if (psize == sizeof(unsigned int))
289  {
290  unsigned int* pdata = (unsigned int*)data;
291 
292  if (pdata[n])
293 
294  ofs << i << '\t' << j << '\t' << k << '\t' << pdata[n] << G4endl;
295 
296  }
297 
298  else if (psize == sizeof(G4double))
299 
300  {
301  G4double* pdata = (G4double*)data;
302  if (pdata[n]) ofs << i << '\t' << j << '\t' << k << '\t' << pdata[n] << G4endl;
303  }
304  }
305 
306  ofs.close();
307  }
308  }
309 }
310 
312 // Store fluence per single ion in multiple files
314 {
315  for (size_t i=0; i < ionStore.size(); i++){
316  StoreMatrix(ionStore[i].name + "_Fluence.out", ionStore[i].fluence, sizeof(unsigned int));
317  }
318 }
319 
321 // Store dose per single ion in multiple files
323 {
324 
325  for (size_t i=0; i < ionStore.size(); i++){
326  StoreMatrix(ionStore[i].name + "_Dose.out", ionStore[i].dose, sizeof(G4double));
327  }
328 }
329 
330 
332 // Store dose into a single file
333 // or in histograms. Please note that this function is called via
334 // messenger commands
335 // defined in the HadrontherapyAnalysisFileMessenger.cc class file
337 {
338 #define width 15L
339  filename = (file=="") ? stdFile:file;
340 
341  // Sort like periodic table
342 
343  std::sort(ionStore.begin(), ionStore.end());
344  G4cout << "Dose is being written to " << filename << G4endl;
345  ofs.open(filename, std::ios::out);
346 
347 
348  if (ofs.is_open())
349  {
350  // Write the voxels index and the list of particles/ions
351 
352  ofs << std::setprecision(6) << std::left <<
353  "i\tj\tk\t";
354 
355 
356  // Total dose
357  ofs << std::setw(width) << "Dose(Gy)";
358 
359 
360  if (secondary)
361  {
362  for (size_t l=0; l < ionStore.size(); l++)
363  {
364  G4String a = (ionStore[l].isPrimary) ? "_1":""; // is it a primary?
365  ofs << std::setw(width) << ionStore[l].name + a <<
366  std::setw(width) << ionStore[l].name + a;
367 
368 
369  }
370  ofs << G4endl;
371  }
372 
373  // Write data
374  for(G4int i = 0; i < numberOfVoxelAlongX; i++)
375  for(G4int j = 0; j < numberOfVoxelAlongY; j++)
376  for(G4int k = 0; k < numberOfVoxelAlongZ; k++)
377  {
378  G4int n = Index(i, j, k);
379 
380  if (matrix[n])
381  {
382  ofs << G4endl;
383  ofs << i << '\t' << j << '\t' << k << '\t';
384 
385  // Total dose
386  ofs << std::setw(width) << (matrix[n]/massOfVoxel)/doseUnit;
387 
388 
389  if (secondary)
390  {
391  for (size_t l=0; l < ionStore.size(); l++)
392  {
393  // Fill ASCII file rows
394  ofs << std::setw(width) << ionStore[l].dose[n]/massOfVoxel/doseUnit <<
395  std::setw(width) << ionStore[l].fluence[n];
396 
397  }
398  }
399  }
400  }
401  ofs.close();
402  }
403 
404 
405 
406 }
409  G4double energyDeposit)
410 {
411  if (matrix)
412  matrix[Index(i,j,k)] += energyDeposit;
413 
414  // Store the energy deposit in the matrix element corresponding
415  // to the phantom voxel
416 }
417 
418 
419 
420