Geant4  10.01
IORTMatrix.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 // This is the *BASIC* version of IORT, a Geant4-based application
27 //
28 // Main Authors: G.Russo(a,b), C.Casarino*(c), G.C. Candiano(c), G.A.P. Cirrone(d), F.Romano(d)
29 // Contributor Authors: S.Guatelli(e)
30 // Past Authors: G.Arnetta(c), S.E.Mazzaglia(d)
31 //
32 // (a) Fondazione Istituto San Raffaele G.Giglio, Cefalù, Italy
33 // (b) IBFM-CNR , Segrate (Milano), Italy
34 // (c) LATO (Laboratorio di Tecnologie Oncologiche), Cefalù, Italy
35 // (d) Laboratori Nazionali del Sud of the INFN, Catania, Italy
36 // (e) University of Wallongong, Australia
37 //
38 // *Corresponding author, email to carlo.casarino@polooncologicocefalu.it
40 
41 #include <fstream>
42 #include <iostream>
43 #include <sstream>
44 #include <iomanip>
45 
46 #include "IORTMatrix.hh"
47 #include "IORTAnalysisManager.hh"
49 
50 #include "globals.hh"
51 #include "G4SystemOfUnits.hh"
52 #include "G4RunManager.hh"
53 #include "G4ParticleGun.hh"
54 
55 // Units definition: CLHEP/Units/SystemOfUnits.h
56 //
59 
60 // Only return a pointer to matrix
62 {
63  return instance;
64 }
65  // This STATIC method delete (!) the old matrix and rewrite a new object returning a pointer to it
66  // TODO A check on the parameters is required!
68 {
69  if (instance) delete instance;
70  instance = new IORTMatrix(voxelX, voxelY, voxelZ, mass);
71  instance -> Initialize();
72  return instance;
73 }
74 IORTMatrix::IORTMatrix(G4int voxelX, G4int voxelY, G4int voxelZ, G4double mass):
75  stdFile("Dose.out"),
76  doseUnit(MeV/g)
77 {
78  // Number of the voxels of the phantom
79  // For Y = Z = 1 the phantom is divided in slices (and not in voxels)
80  // orthogonal to the beam axis
81  numberOfVoxelAlongX = voxelX;
82  numberOfVoxelAlongY = voxelY;
83  numberOfVoxelAlongZ = voxelZ;
84  massOfVoxel = mass;
85  // Create the dose matrix
87  if (matrix)
88  {
89  G4cout << "IORTMatrix: Memory space to store physical dose into " <<
91  " voxels has been allocated " << G4endl;
92  }
93  else G4Exception("IORTMatrix::IORTMatrix()", "IORT0005", FatalException, "Error: can't allocate memory to store physical dose!");
94  // Hit voxel (TrackID) marker
95  // This array mark the status of voxel, if a hit occur, with the trackID of the particle
96  // Must be initialized
98  ClearHitTrack();
99 }
100 
103 {
104  delete[] matrix;
105  delete[] hitTrack;
106  // free fluences/dose data memory
107  Clear();
108 }
109 
112 {
113  for (size_t i=0; i<ionStore.size(); i++)
114  {
115  delete[] ionStore[i].dose;
116  delete[] ionStore[i].fluence;
117  }
118  ionStore.clear();
119 }
120 
122 // Initialise the elements of the matrix to zero
124 {
125  // Clear ions store
126  Clear();
127  // Clear dose
129  {
130  matrix[i] = 0;
131  }
132 }
135  // Print generated nuclides list
137 {
138  for (size_t i=0; i<ionStore.size(); i++)
139  {
140  G4cout << ionStore[i].name << G4endl;
141  }
142 }
144  // Clear Hit voxel (TrackID) markers
146 {
148 }
149  // Return Hit status
151 {
152  return &(hitTrack[Index(i,j,k)]);
153 }
155 // Dose methods...
156 // Fill DOSE/fluence matrix for secondary particles:
157 // If fluence parameter is true (default value is FALSE) then fluence at voxel (i, j, k) is increased.
158 // The energyDeposit parameter fill the dose matrix for voxel (i,j,k)
160 
162  G4ParticleDefinition* particleDef,
163  G4int i, G4int j, G4int k,
164  G4double energyDeposit,
165  G4bool fluence)
166 {
167  if ( (energyDeposit <=0. && !fluence) || !secondary) return false;
168  // Get Particle Data Group particle ID
169  G4int PDGencoding = particleDef -> GetPDGEncoding();
170  PDGencoding -= PDGencoding%10;
171 
172  // Search for already allocated data...
173  for (size_t l=0; l < ionStore.size(); l++)
174  {
175  if (ionStore[l].PDGencoding == PDGencoding )
176  { // Is it a primary or a secondary particle?
177  if ( ((trackID == 1) && (ionStore[l].isPrimary)) || ((trackID !=1) && (!ionStore[l].isPrimary)))
178  {
179  if (energyDeposit > 0.) ionStore[l].dose[Index(i, j, k)] += energyDeposit/massOfVoxel;
180 
181  // Fill a matrix per each ion with the fluence
182  if (fluence) ionStore[l].fluence[Index(i, j, k)]++;
183  return true;
184  }
185  }
186  }
187 
188  G4int Z = particleDef-> GetAtomicNumber();
189  G4int A = particleDef-> GetAtomicMass();
190 
191  G4String fullName = particleDef -> GetParticleName();
192  G4String name = fullName.substr (0, fullName.find("[") ); // cut excitation energy
193  // Let's put a new particle in our store...
194  ion newIon =
195  {
196  (trackID == 1) ? true:false,
197  PDGencoding,
198  name,
199  name.length(),
200  Z,
201  A,
204  };
205  // Initialize data
206  if (newIon.dose && newIon.fluence)
207  {
209  {
210  newIon.dose[q] = 0.;
211  newIon.fluence[q] = 0;
212  }
213  if (energyDeposit > 0.) newIon.dose[Index(i, j, k)] += energyDeposit/massOfVoxel;
214  if (fluence) newIon.fluence[Index(i, j, k)]++;
215 
216  ionStore.push_back(newIon);
217 
218  // TODO Put some verbosity check
219  /*
220  G4cout << "Memory space to store the DOSE/FLUENCE into " <<
221  numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ <<
222  " voxels has been allocated for the nuclide " << newIon.name <<
223  " (Z = " << Z << ", A = " << A << ")" << G4endl ;
224  */
225  return true;
226  }
227  else // XXX Out of memory! XXX
228  {
229  return false;
230  }
231 
232 }
233 
236  // Methods to store data to filenames...
239  //
240  // General method to store matrix data to filename
241 void IORTMatrix::StoreMatrix(G4String file, void* data, size_t psize)
242 {
243  if (data)
244  {
245  ofs.open(file, std::ios::out);
246  if (ofs.is_open())
247  {
248  for(G4int i = 0; i < numberOfVoxelAlongX; i++)
249  for(G4int j = 0; j < numberOfVoxelAlongY; j++)
250  for(G4int k = 0; k < numberOfVoxelAlongZ; k++)
251  {
252  G4int n = Index(i, j, k);
253  // Check for data type: u_int, G4double, XXX
254  if (psize == sizeof(unsigned int))
255  {
256  unsigned int* pdata = (unsigned int*)data;
257  if (pdata[n]) ofs << i << '\t' << j << '\t' <<
258  k << '\t' << pdata[n] << G4endl;
259  }
260  else if (psize == sizeof(G4double))
261  {
262  G4double* pdata = (G4double*)data;
263  if (pdata[n]) ofs << i << '\t' << j << '\t' <<
264  k << '\t' << pdata[n] << G4endl;
265  }
266  }
267  ofs.close();
268  }
269  }
270 }
271 
272  // Store fluence per single ion in multiple files
274 {
275  for (size_t i=0; i < ionStore.size(); i++){
276  StoreMatrix(ionStore[i].name + "_Fluence.out", ionStore[i].fluence, sizeof(unsigned int));
277  }
278 }
279  // Store dose per single ion in multiple files
281 {
282 
283  for (size_t i=0; i < ionStore.size(); i++){
284  StoreMatrix(ionStore[i].name + "_Dose.out", ionStore[i].dose, sizeof(G4double));
285  }
286 }
288  // Store dose for all ions into a single file and into ntuples.
289  // Please note that this function is called via messenger commands
290  // defined in the IORTAnalysisFileMessenger.cc class file
292 {
293 #define width 15L
294  filename = (file=="") ? stdFile:file;
295  // Sort like periodic table
296  std::sort(ionStore.begin(), ionStore.end());
297  G4cout << "Dose is being written to " << filename << G4endl;
298  ofs.open(filename, std::ios::out);
299  if (ofs.is_open())
300  {
301  // Write the voxels index and the list of particles/ions
302  ofs << std::setprecision(6) << std::left <<
303  "i\tj\tk\t";
304  // Total dose
305  ofs << std::setw(width) << "Dose(MeV/g)";
306  if (secondary)
307  {
308  for (size_t l=0; l < ionStore.size(); l++)
309  {
310  G4String a = (ionStore[l].isPrimary) ? "_1":""; // is it a primary?
311  ofs << std::setw(width) << ionStore[l].name + a <<
312  std::setw(width) << ionStore[l].name + a;
313  }
314  ofs << G4endl;
315 
316  /*
317  * PDGencondig
318  */
319  /*
320  ofs << std::setprecision(6) << std::left <<
321  "0\t0\t0\t";
322 
323  // Total dose
324  ofs << std::setw(width) << '0';
325  for (size_t l=0; l < ionStore.size(); l++)
326  {
327  ofs << std::setw(width) << ionStore[l].PDGencoding <<
328  std::setw(width) << ionStore[l].PDGencoding;
329  }
330  ofs << G4endl;
331  */
332  }
333  // Write data
334  for(G4int i = 0; i < numberOfVoxelAlongX; i++)
335  for(G4int j = 0; j < numberOfVoxelAlongY; j++)
336  for(G4int k = 0; k < numberOfVoxelAlongZ; k++)
337  {
338  G4int n = Index(i, j, k);
339  // Write only not identically null data lines
340  if (matrix[n])
341  {
342  ofs << G4endl;
343  ofs << i << '\t' << j << '\t' << k << '\t';
344  // Total dose
345  ofs << std::setw(width) << matrix[n]/massOfVoxel/doseUnit;
346  if (secondary)
347  {
348  for (size_t l=0; l < ionStore.size(); l++)
349  {
350  // Fill ASCII file rows
351  ofs << std::setw(width) << ionStore[l].dose[n]/massOfVoxel/doseUnit <<
352  std::setw(width) << ionStore[l].fluence[n];
353  }
354  }
355  }
356  }
357  ofs.close();
358  }
359 }
361 
362 #ifdef G4ANALYSIS_USE_ROOT
363 void IORTMatrix::StoreDoseFluenceRoot()
364 {
366  if (analysis -> IsTheTFile())
367  {
368  for(G4int i = 0; i < numberOfVoxelAlongX; i++)
369  for(G4int j = 0; j < numberOfVoxelAlongY; j++)
370  for(G4int k = 0; k < numberOfVoxelAlongZ; k++)
371  {
372  G4int n = Index(i, j, k);
373  for (size_t l=0; l < ionStore.size(); l++)
374 
375  {
376  // Do the same work for .root file: fill dose/fluence ntuple
377  analysis -> FillVoxelFragmentTuple( i, j, k,
378  ionStore[l].A,
379  ionStore[l].Z,
380  ionStore[l].dose[n]/massOfVoxel/doseUnit,
381  ionStore[l].fluence[n] );
382 
383 
384  }
385  }
386  }
387 }
388 #endif
389 
391  G4double energyDeposit)
392 {
393  if (matrix)
394  matrix[Index(i,j,k)] += energyDeposit;
395 
396  // Store the energy deposit in the matrix element corresponding
397  // to the phantom voxel
398 }
400 {
401  // Store the information of the matrix in a ntuple and in
402  // a 1D Histogram
403 
404 /*
406  G4int k;
407  G4int j;
408  G4int i;
409 
410  if (matrix)
411  { // AGGIUNTO
412  std::ofstream ofs; // AGGIUNTO
413 
414 ofs.open("PDD9.9Mev_coll60_0gradi_s500_Sp1_6gradi_step0.01_setcuts0.01.out"); // AGGIUNTO
415 
416  for(G4int l = 0; l < numberOfVoxelAlongZ; l++) // was "numberVoxelZ" and so in the other directions
417  {
418  k = l;
419 
420  for(G4int m = 0; m < numberOfVoxelAlongY; m++)
421  {
422  j = m * numberOfVoxelAlongZ + k;
423 
424  for(G4int n = 0; n < numberOfVoxelAlongX; n++)
425  {
426  i = n* numberOfVoxelAlongZ * numberOfVoxelAlongY + j;
427  if(matrix[i] != 0)
428  {
429  ofs<< n <<'\t'<< m <<'\t'<< // AGGIUNTO
430  k<<'\t'<<matrix[i]<<G4endl; // AGGIUNTO
431 
432 
433  }
434  }
435  }
436  }
437  ofs.close(); // AGGIUNTO
438  }
440 */
441 
442  // Convert energy deposited to dose.
443  // Store the information of the matrix in a ntuple and in
444  // a 1D Histogram
445 /*
446  IORTAnalysisManager* analysis = IORTAnalysisManager::GetInstance();
447  if (matrix)
448  {
449  for(G4int i = 0; i < numberOfVoxelAlongX; i++)
450  for(G4int j = 0; j < numberOfVoxelAlongY; j++)
451  for(G4int k = 0; k < numberOfVoxelAlongZ; k++)
452  {
453  G4int n = Index(i,j,k);
454 
455 #ifdef G4ANALYSIS_USE_ROOT
456  if (analysis -> IsTheTFile() )
457  {
458  analysis -> FillEnergyDeposit(i, j, k, matrix[n]/massOfVoxel/doseUnit);
459  analysis -> BraggPeak(i, matrix[n]/massOfVoxel/doseUnit);
460  }
461 #endif
462 
463  }
464  }
465 */
466 }
467 
468 
469 
470 
static const double MeV
Definition: G4SIunits.hh:193
void Clear()
Definition: IORTMatrix.cc:111
G4int numberOfVoxelAlongY
Definition: IORTMatrix.hh:136
G4double * dose
G4String name
Definition: TRTMaterials.hh:40
void StoreMatrix(G4String file, void *data, size_t psize)
Definition: IORTMatrix.cc:241
static G4bool secondary
Definition: IORTMatrix.hh:83
G4double a
Definition: TRTMaterials.hh:39
G4int * GetHitTrack(G4int i, G4int j, G4int k)
Definition: IORTMatrix.cc:150
G4double massOfVoxel
Definition: IORTMatrix.hh:138
A class for connecting the simulation to an analysis package.
int G4int
Definition: G4Types.hh:78
G4double doseUnit
Definition: IORTMatrix.hh:148
void TotalEnergyDeposit()
Definition: IORTMatrix.cc:399
G4String stdFile
Definition: IORTMatrix.hh:142
G4GLOB_DLL std::ostream G4cout
G4double * matrix
Definition: IORTMatrix.hh:140
unsigned int * fluence
bool G4bool
Definition: G4Types.hh:79
void StoreDoseData()
Definition: IORTMatrix.cc:280
void StoreFluenceData()
Definition: IORTMatrix.cc:273
G4int numberOfVoxelAlongX
Definition: IORTMatrix.hh:135
void ClearHitTrack()
Definition: IORTMatrix.cc:145
const G4int n
static const G4double A[nN]
static IORTAnalysisManager * GetInstance()
Get the pointer to the analysis manager.
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4bool Fill(G4int, G4ParticleDefinition *particleDef, G4int i, G4int j, G4int k, G4double energyDeposit, G4bool fluence=false)
Definition: IORTMatrix.cc:161
static IORTMatrix * instance
Definition: IORTMatrix.hh:134
static const double g
Definition: G4SIunits.hh:162
IORTMatrix(G4int numberOfVoxelAlongX, G4int numberOfVoxelAlongY, G4int numberOfVoxelAlongZ, G4double massOfVoxel)
Definition: IORTMatrix.cc:74
void PrintNuclides()
Definition: IORTMatrix.cc:136
G4int * hitTrack
Definition: IORTMatrix.hh:141
void Initialize()
Definition: IORTMatrix.cc:123
#define G4endl
Definition: G4ios.hh:61
static IORTMatrix * GetInstance()
Definition: IORTMatrix.cc:61
std::vector< ion > ionStore
Definition: IORTMatrix.hh:146
double G4double
Definition: G4Types.hh:76
G4String filename
Definition: IORTMatrix.hh:142
G4int numberOfVoxelAlongZ
Definition: IORTMatrix.hh:137
std::ofstream ofs
Definition: IORTMatrix.hh:143
#define width
void StoreDoseFluenceAscii(G4String filename="")
Definition: IORTMatrix.cc:291
G4int Index(G4int i, G4int j, G4int k)
Definition: IORTMatrix.hh:122