Geant4  10.02.p01
HadrontherapyMatrix.cc
Go to the documentation of this file.
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"
37 #include "globals.hh"
38 #include "G4SystemOfUnits.hh"
39 #include "G4RunManager.hh"
40 #include "G4ParticleGun.hh"
41 
42 // Units definition: CLHEP/Units/SystemOfUnits.h
43 //
46 
47 // Only return a pointer to matrix
49 {
50  return instance;
51 }
52  // This STATIC method delete (!) the old matrix and rewrite a new object returning a pointer to it
53  // TODO A check on the parameters is required!
55 {
56  if (instance) delete instance;
57  instance = new HadrontherapyMatrix(voxelX, voxelY, voxelZ, mass);
58  instance -> Initialize();
59  return instance;
60 }
62  stdFile("Dose.out"),
63  doseUnit(gray)
64 {
65  // Number of the voxels of the phantom
66  // For Y = Z = 1 the phantom is divided in slices (and not in voxels)
67  // orthogonal to the beam axis
68  numberOfVoxelAlongX = voxelX;
69  numberOfVoxelAlongY = voxelY;
70  numberOfVoxelAlongZ = voxelZ;
71  massOfVoxel = mass;
72  // Create the dose matrix
74  if (matrix)
75  {
76  G4cout << "HadrontherapyMatrix: Memory space to store physical dose into " <<
78  " voxels has been allocated " << G4endl;
79  }
80  else G4Exception("HadrontherapyMatrix::HadrontherapyMatrix()", "Hadrontherapy0005", FatalException, "Can't allocate memory to store physical dose!");
81  // Hit voxel (TrackID) marker
82  // This array mark the status of voxel, if a hit occur, with the trackID of the particle
83  // Must be initialized
85  ClearHitTrack();
86 }
87 
90 {
91  delete[] matrix;
92  delete[] hitTrack;
93  // free fluences/dose data memory
94  Clear();
95 }
96 
99 {
100  for (size_t i=0; i<ionStore.size(); i++)
101  {
102  delete[] ionStore[i].dose;
103  delete[] ionStore[i].fluence;
104  }
105  ionStore.clear();
106 }
107 
109 // Initialise the elements of the matrix to zero
111 {
112  // Clear ions store
113  Clear();
114  // Clear dose
116  {
117  matrix[i] = 0;
118  }
119 }
122  // Print generated nuclides list
124 {
125  for (size_t i=0; i<ionStore.size(); i++)
126  {
127  G4cout << ionStore[i].name << G4endl;
128  }
129 }
131  // Clear Hit voxel (TrackID) markers
133 {
135 }
136  // Return Hit status
138 {
139  return &(hitTrack[Index(i,j,k)]);
140 }
142 // Dose methods...
143 // Fill DOSE/fluence matrix for secondary particles:
144 // If fluence parameter is true (default value is FALSE) then fluence at voxel (i, j, k) is increased.
145 // The energyDeposit parameter fill the dose matrix for voxel (i,j,k)
147 
149  G4ParticleDefinition* particleDef,
150  G4int i, G4int j, G4int k,
151  G4double energyDeposit,
152  G4bool fluence)
153 {
154  if ( (energyDeposit <=0. && !fluence) || !secondary) return false;
155  // Get Particle Data Group particle ID
156  G4int PDGencoding = particleDef -> GetPDGEncoding();
157  PDGencoding -= PDGencoding%10;
158 
159  // Search for already allocated data...
160  for (size_t l=0; l < ionStore.size(); l++)
161  {
162  if (ionStore[l].PDGencoding == PDGencoding )
163  { // Is it a primary or a secondary particle?
164  if ( (trackID ==1 && ionStore[l].isPrimary) || (trackID !=1 && !ionStore[l].isPrimary))
165  {
166  if (energyDeposit > 0.) ionStore[l].dose[Index(i, j, k)] += energyDeposit;
167 
168  // Fill a matrix per each ion with the fluence
169  if (fluence) ionStore[l].fluence[Index(i, j, k)]++;
170  return true;
171  }
172  }
173  }
174 
175  G4int Z = particleDef-> GetAtomicNumber();
176  G4int A = particleDef-> GetAtomicMass();
177 
178  G4String fullName = particleDef -> GetParticleName();
179  G4String name = fullName.substr (0, fullName.find("[") ); // cut excitation energy
180  // Let's put a new particle in our store...
181  ion newIon =
182  {
183  (trackID == 1) ? true:false,
184  PDGencoding,
185  name,
186  name.length(),
187  Z,
188  A,
191  };
192  // Initialize data
193  if (newIon.dose && newIon.fluence)
194  {
196  {
197  newIon.dose[q] = 0.;
198  newIon.fluence[q] = 0;
199  }
200  if (energyDeposit > 0.) newIon.dose[Index(i, j, k)] += energyDeposit;
201  if (fluence) newIon.fluence[Index(i, j, k)]++;
202 
203  ionStore.push_back(newIon);
204 
205  // TODO Put some verbosity check
206  /*
207  G4cout << "Memory space to store the DOSE/FLUENCE into " <<
208  numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ <<
209  " voxels has been allocated for the nuclide " << newIon.name <<
210  " (Z = " << Z << ", A = " << A << ")" << G4endl ;
211  */
212  return true;
213  }
214  else // XXX Out of memory! XXX
215  {
216  return false;
217  }
218 
219 }
220 
223  // Methods to store data to filenames...
226  //
227  // General method to store matrix data to filename
228 void HadrontherapyMatrix::StoreMatrix(G4String file, void* data, size_t psize)
229 {
230  if (data)
231  {
232  ofs.open(file, std::ios::out);
233  if (ofs.is_open())
234  {
235  for(G4int i = 0; i < numberOfVoxelAlongX; i++)
236  for(G4int j = 0; j < numberOfVoxelAlongY; j++)
237  for(G4int k = 0; k < numberOfVoxelAlongZ; k++)
238  {
239  G4int n = Index(i, j, k);
240  // Check for data type: u_int, G4double, XXX
241  if (psize == sizeof(unsigned int))
242  {
243  unsigned int* pdata = (unsigned int*)data;
244  if (pdata[n]) ofs << i << '\t' << j << '\t' <<
245  k << '\t' << pdata[n] << G4endl;
246  }
247  else if (psize == sizeof(G4double))
248  {
249  G4double* pdata = (G4double*)data;
250  if (pdata[n]) ofs << i << '\t' << j << '\t' <<
251  k << '\t' << pdata[n] << G4endl;
252  }
253  }
254  ofs.close();
255  }
256  }
257 }
258 
259  // Store fluence per single ion in multiple files
261 {
262  for (size_t i=0; i < ionStore.size(); i++){
263  StoreMatrix(ionStore[i].name + "_Fluence.out", ionStore[i].fluence, sizeof(unsigned int));
264  }
265 }
266  // Store dose per single ion in multiple files
268 {
269 
270  for (size_t i=0; i < ionStore.size(); i++){
271  StoreMatrix(ionStore[i].name + "_Dose.out", ionStore[i].dose, sizeof(G4double));
272  }
273 }
275  // Store dose for all ions into a single file and into ntuples.
276  // Please note that this function is called via messenger commands
277  // defined in the HadrontherapyAnalysisFileMessenger.cc class file
279 {
280 #define width 15L
281  filename = (file=="") ? stdFile:file;
282  // Sort like periodic table
283  std::sort(ionStore.begin(), ionStore.end());
284  G4cout << "Dose is being written to " << filename << G4endl;
285  ofs.open(filename, std::ios::out);
286  if (ofs.is_open())
287  {
288  // Write the voxels index and the list of particles/ions
289  ofs << std::setprecision(6) << std::left <<
290  "i\tj\tk\t";
291  // Total dose
292  ofs << std::setw(width) << "Dose(Gy)";
293  if (secondary)
294  {
295  for (size_t l=0; l < ionStore.size(); l++)
296  {
297  G4String a = (ionStore[l].isPrimary) ? "_1":""; // is it a primary?
298  ofs << std::setw(width) << ionStore[l].name + a <<
299  std::setw(width) << ionStore[l].name + a;
300  }
301  ofs << G4endl;
302 
303  /*
304  * PDGencondig
305  */
306  /*
307  ofs << std::setprecision(6) << std::left <<
308  "0\t0\t0\t";
309 
310  // Total dose
311  ofs << std::setw(width) << '0';
312  for (size_t l=0; l < ionStore.size(); l++)
313  {
314  ofs << std::setw(width) << ionStore[l].PDGencoding <<
315  std::setw(width) << ionStore[l].PDGencoding;
316  }
317  ofs << G4endl;
318  */
319  }
320  // Write data
321  for(G4int i = 0; i < numberOfVoxelAlongX; i++)
322  for(G4int j = 0; j < numberOfVoxelAlongY; j++)
323  for(G4int k = 0; k < numberOfVoxelAlongZ; k++)
324  {
325  G4int n = Index(i, j, k);
326  // Write only not identically null data lines
327  if (matrix[n])
328  {
329  ofs << G4endl;
330  ofs << i << '\t' << j << '\t' << k << '\t';
331  // Total dose
332  ofs << std::setw(width) << (matrix[n]/massOfVoxel)/doseUnit;
333  if (secondary)
334  {
335  for (size_t l=0; l < ionStore.size(); l++)
336  {
337  // Fill ASCII file rows
338  ofs << std::setw(width) << ionStore[l].dose[n]/massOfVoxel/doseUnit <<
339  std::setw(width) << ionStore[l].fluence[n];
340  }
341  }
342  }
343  }
344  ofs.close();
345  }
346 }
348 
349 #ifdef G4ANALYSIS_USE_ROOT
350 void HadrontherapyMatrix::StoreDoseFluenceRoot()
351 {
353  if (analysis -> IsTheTFile())
354  {
355  for(G4int i = 0; i < numberOfVoxelAlongX; i++)
356  for(G4int j = 0; j < numberOfVoxelAlongY; j++)
357  for(G4int k = 0; k < numberOfVoxelAlongZ; k++)
358  {
359  G4int n = Index(i, j, k);
360  for (size_t l=0; l < ionStore.size(); l++)
361 
362  {
363  // Do the same work for .root file: fill dose/fluence ntuple
364  analysis -> FillVoxelFragmentTuple( i, j, k,
365  ionStore[l].A,
366  ionStore[l].Z,
367  ionStore[l].dose[n]/massOfVoxel/doseUnit,
368  ionStore[l].fluence[n] );
369 
370 
371  }
372  }
373  }
374 }
375 #endif
376 
378  G4double energyDeposit)
379 {
380  if (matrix)
381  matrix[Index(i,j,k)] += energyDeposit;
382 
383  // Store the energy deposit in the matrix element corresponding
384  // to the phantom voxel
385 }
387 {
388  // Convert energy deposited to dose.
389  // Store the information of the matrix in a ntuple and in
390  // a 1D Histogram
391 #ifdef G4ANALYSIS_USE_ROOT
393 #endif
394 
395  if (matrix)
396  {
397  for(G4int i = 0; i < numberOfVoxelAlongX; i++)
398  for(G4int j = 0; j < numberOfVoxelAlongY; j++)
399  for(G4int k = 0; k < numberOfVoxelAlongZ; k++)
400  {
401 #ifdef G4ANALYSIS_USE_ROOT
402  G4int n = Index(i,j,k);
403  if (analysis -> IsTheTFile() )
404  {
405  analysis -> FillEnergyDeposit(i, j, k, matrix[n]/massOfVoxel/doseUnit);
406  analysis -> BraggPeak(i, matrix[n]/massOfVoxel/doseUnit);
407  }
408 #endif
409  }
410  }
411 }
412 
static HadrontherapyAnalysisManager * GetInstance()
Get the pointer to the analysis manager.
G4double * dose
G4int Index(G4int i, G4int j, G4int k)
G4String name
Definition: TRTMaterials.hh:40
#define width
void StoreDoseFluenceAscii(G4String filename="")
G4double a
Definition: TRTMaterials.hh:39
int G4int
Definition: G4Types.hh:78
static HadrontherapyMatrix * GetInstance()
void StoreMatrix(G4String file, void *data, size_t psize)
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
unsigned int * fluence
bool G4bool
Definition: G4Types.hh:79
G4int * GetHitTrack(G4int i, G4int j, G4int k)
G4bool Fill(G4int, G4ParticleDefinition *particleDef, G4int i, G4int j, G4int k, G4double energyDeposit, G4bool fluence=false)
const G4int n
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static const double gray
Definition: G4SIunits.hh:306
HadrontherapyMatrix(G4int numberOfVoxelAlongX, G4int numberOfVoxelAlongY, G4int numberOfVoxelAlongZ, G4double massOfVoxel)
A class for connecting the simulation to an analysis package.
static HadrontherapyMatrix * instance
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
std::vector< ion > ionStore