Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
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"
41 
42 
44 #include "G4SystemOfUnits.hh"
45 #include <time.h>
46 
49 
50 
51 {
53 }
55 {
56  delete fMess;
57 }
58 
60 
62  return instance;
63 }
64 
65 
66 
69 
70 
71 // Only return a pointer to matrix
73 {
74  return instance;
75 }
76 // This STATIC method delete (!) the old matrix and rewrite a new object returning a pointer to it
77 // TODO A check on the parameters is required!
79 {
80  if (instance) delete instance;
81  instance = new HadrontherapyMatrix(voxelX, voxelY, voxelZ, mass);
82  instance -> Initialize();
83  return instance;
84 }
85 
86 
88 stdFile("Dose.out"),
89 doseUnit(gray)
90 
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 
158 
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)
190 
192  G4ParticleDefinition* particleDef,
193  G4int i, G4int j, G4int k,
194  G4double energyDeposit,
195  G4bool fluence)
196 {
197 
198  if ( (energyDeposit <=0. && !fluence) || !secondary) return false;
199 
200  // Get Particle Data Group particle ID
201  G4int PDGencoding = particleDef -> GetPDGEncoding();
202  PDGencoding -= PDGencoding%10;
203 
204  // Search for already allocated data...
205  for (size_t l=0; l < ionStore.size(); l++)
206  {
207  if (ionStore[l].PDGencoding == PDGencoding )
208  { // Is it a primary or a secondary particle?
209 
210  if ( (trackID ==1 && ionStore[l].isPrimary) || (trackID !=1 && !ionStore[l].isPrimary))
211  {
212  if (energyDeposit > 0.)
213 
214  ionStore[l].dose[Index(i, j, k)] += energyDeposit;
215 
216  // Fill a matrix per each ion with the fluence
217 
218  if (fluence) ionStore[l].fluence[Index(i, j, k)]++;
219  return true;
220  }
221  }
222  }
223 
224  G4int Z = particleDef-> GetAtomicNumber();
225  G4int A = particleDef-> GetAtomicMass();
226  G4String fullName = particleDef -> GetParticleName();
227  G4String name = fullName.substr (0, fullName.find("[") ); // cut excitation energy
228 
229  // Let's put a new particle in our store...
230 
231 
232  ion newIon =
233  {
234  (trackID == 1) ? true:false,
235  PDGencoding,
236  name,
237  name.length(),
238  Z,
239  A,
242  };
243 
244 
245  // Initialize data
246  if (newIon.dose && newIon.fluence)
247  {
249  {
250  newIon.dose[q] = 0.;
251  newIon.fluence[q] = 0;
252  }
253 
254  if (energyDeposit > 0.) newIon.dose[Index(i, j, k)] += energyDeposit;
255  if (fluence) newIon.fluence[Index(i, j, k)]++;
256 
257  ionStore.push_back(newIon);
258  return true;
259  }
260 
261  else // XXX Out of memory! XXX
262 
263  {
264  return false;
265  }
266 
267 }
268 
271 // Methods to store data to filenames...
274 //
275 // General method to store matrix data to filename
276 
277 
279 {
280  if (data)
281  {
282  ofs.open(file, std::ios::out);
283  if (ofs.is_open())
284  {
285  for(G4int i = 0; i < numberOfVoxelAlongX; i++)
286  for(G4int j = 0; j < numberOfVoxelAlongY; j++)
287  for(G4int k = 0; k < numberOfVoxelAlongZ; k++)
288  {
289  G4int n = Index(i, j, k);
290 
291  if (psize == sizeof(unsigned int))
292  {
293  unsigned int* pdata = (unsigned int*)data;
294 
295  if (pdata[n])
296 
297  ofs << i << '\t' << j << '\t' << k << '\t' << pdata[n] << G4endl;
298 
299  }
300 
301  else if (psize == sizeof(G4double))
302 
303  {
304  G4double* pdata = (G4double*)data;
305  if (pdata[n]) ofs << i << '\t' << j << '\t' << k << '\t' << pdata[n] << G4endl;
306  }
307  }
308 
309  ofs.close();
310  }
311  }
312 }
313 
314 // Store fluence per single ion in multiple files
316 {
317  for (size_t i=0; i < ionStore.size(); i++){
318  StoreMatrix(ionStore[i].name + "_Fluence.out", ionStore[i].fluence, sizeof(unsigned int));
319  }
320 }
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 for all ions into a single file and into ntuples.
333 // Please note that this function is called via messenger commands
334 // defined in the HadrontherapyAnalysisFileMessenger.cc class file
335 
336 
338 {
339 
340 
341 #define width 15L
342  filename = (file=="") ? stdFile:file;
343 
344  // Sort like periodic table
345 
346  std::sort(ionStore.begin(), ionStore.end());
347  G4cout << "Dose is being written to " << filename << G4endl;
348  ofs.open(filename, std::ios::out);
349 
350 
351  if (ofs.is_open())
352  {
353  // Write the voxels index and the list of particles/ions
354 
355  ofs << std::setprecision(6) << std::left <<
356  "i\tj\tk\t";
357 
358 
359  // Total dose
360  ofs << std::setw(width) << "Dose(Gy)";
361 
362 
363  if (secondary)
364  {
365  for (size_t l=0; l < ionStore.size(); l++)
366  {
367  G4String a = (ionStore[l].isPrimary) ? "_1":""; // is it a primary?
368  ofs << std::setw(width) << ionStore[l].name + a <<
369  std::setw(width) << ionStore[l].name + a;
370 
371 
372  }
373  ofs << G4endl;
374  }
375 
376 
377  G4AnalysisManager* Analysis = G4AnalysisManager::Instance();
378 
379  Analysis ->SetVerboseLevel(1);
380  Analysis ->SetFirstHistoId(1);
381  Analysis ->SetFirstNtupleId(1);
382  Analysis ->OpenFile("Dose");
383 
384 
385  Analysis ->CreateNtuple("coordinate", "dose");
386 
387  Analysis ->CreateNtupleIColumn("i");//1
388  Analysis ->CreateNtupleIColumn("j");//2
389  Analysis ->CreateNtupleIColumn("k");//3
390  Analysis ->CreateNtupleDColumn("totaldose");//6
391  Analysis ->CreateNtupleIColumn("A");//4
392  Analysis ->CreateNtupleIColumn("Z");//5
393  Analysis ->CreateNtupleDColumn("Iondose");//6
394  Analysis ->CreateNtupleDColumn("fluence");//7
395  Analysis ->FinishNtuple();
396 
397 
398 
399  // Write data
400  for(G4int i = 0; i < numberOfVoxelAlongX; i++)
401  for(G4int j = 0; j < numberOfVoxelAlongY; j++)
402  for(G4int k = 0; k < numberOfVoxelAlongZ; k++)
403  {
404  G4int n = Index(i, j, k);
405  // Write only not identically null data lines
406 
407 
408  Analysis->FillNtupleIColumn(1,0, i);
409  Analysis->FillNtupleIColumn(1,1, j);
410  Analysis->FillNtupleIColumn(1,2, k);
411  if (matrix[n])
412  {
413  ofs << G4endl;
414  ofs << i << '\t' << j << '\t' << k << '\t';
415  // Total dose
416  ofs << std::setw(width) << (matrix[n]/massOfVoxel)/doseUnit;
417 
418 
419  Analysis->FillNtupleDColumn(1,3, (matrix[n]/massOfVoxel)/doseUnit);
420  if (secondary)
421  {
422  for (size_t l=0; l < ionStore.size(); l++)
423  {
424  // Fill ASCII file rows
425  ofs << std::setw(width) << ionStore[l].dose[n]/massOfVoxel/doseUnit <<
426  std::setw(width) << ionStore[l].fluence[n];
427 
428 
429  Analysis->FillNtupleIColumn(1,4, ionStore[l].A);
430  Analysis->FillNtupleIColumn(1,5, ionStore[l].Z);
431 
432  Analysis->FillNtupleDColumn(1,6, ionStore[l].dose[n]/massOfVoxel/doseUnit);
433  Analysis->FillNtupleDColumn(1,7, ionStore[l].fluence[n]);
434  Analysis->AddNtupleRow(1);
435 
436 
437 
438  }
439  }
440  }
441  }
442  ofs.close();
443 
444  Analysis->Write();
445  Analysis->CloseFile();
446  }
447 
448 
449 
450 }
452  G4double energyDeposit)
453 {
454  if (matrix)
455  matrix[Index(i,j,k)] += energyDeposit;
456 
457  // Store the energy deposit in the matrix element corresponding
458  // to the phantom voxel
459 }
460 
461 
462 
463 
const XML_Char * name
Definition: expat.h:151
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
static HadrontherapyAnalysisManager * instance
Float_t dose
G4double * dose
static HadrontherapyMatrix * GetInstance()
#define G4endl
Definition: G4ios.hh:61
G4bool SetFirstHistoId(G4int firstId)
G4bool OpenFile(const G4String &fileName="")
G4int Index(G4int i, G4int j, G4int k)
static constexpr double gray
Definition: G4SIunits.hh:309
G4bool FillNtupleIColumn(G4int id, G4int value)
G4int CreateNtupleIColumn(const G4String &name)
G4bool FillNtupleDColumn(G4int id, G4double value)
std::vector< ion > ionStore
const XML_Char const XML_Char * data
Definition: expat.h:268
Float_t Z
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
G4bool Fill(G4int, G4ParticleDefinition *particleDef, G4int i, G4int j, G4int k, G4double energyDeposit, G4bool fluence=false)
unsigned int * fluence
void StoreDoseFluenceAscii(G4String filename="")
double A(double temperature)
G4int CreateNtupleDColumn(const G4String &name)
HadrontherapyAnalysisFileMessenger * fMess
void SetVerboseLevel(G4int verboseLevel)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
#define width
int G4int
Definition: G4Types.hh:78
TFile * file
static HadrontherapyMatrix * instance
G4GLOB_DLL std::ostream G4cout
G4int CreateNtuple(const G4String &name, const G4String &title)
Char_t n[5]
G4int * GetHitTrack(G4int i, G4int j, G4int k)
static HadrontherapyAnalysisManager * GetInstance()
void StoreMatrix(G4String file, void *data, size_t psize)
HadrontherapyMatrix(G4int numberOfVoxelAlongX, G4int numberOfVoxelAlongY, G4int numberOfVoxelAlongZ, G4double massOfVoxel)
G4bool SetFirstNtupleId(G4int firstId)