Geant4  10.01.p01
HadrontherapyLet.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 // ********************************************************************
26 #include "HadrontherapyLet.hh"
30 #include "HadrontherapyMatrix.hh"
31 #include "G4RunManager.hh"
32 #include "G4SystemOfUnits.hh"
33 #include <cmath>
34 
37 
39 {
40  if (instance) delete instance;
41  instance = new HadrontherapyLet(pDet);
42  return instance;
43 }
44 
46 {
47  return instance;
48 }
49 
51  :filename("Let.out"),matrix(0) // Default output filename
52 {
53 
55 
56  if (!matrix)
57  G4Exception("HadrontherapyLet::HadrontherapyLet",
58  "Hadrontherapy0005", FatalException,
59  "HadrontherapyMatrix not found. Firstly create an instance of it.");
60 
61  nVoxels = matrix -> GetNvoxel();
62 
63  numberOfVoxelAlongX = matrix -> GetNumberOfVoxelAlongX();
64  numberOfVoxelAlongY = matrix -> GetNumberOfVoxelAlongY();
65  numberOfVoxelAlongZ = matrix -> GetNumberOfVoxelAlongZ();
66 
68  pPGA = (HadrontherapyPrimaryGeneratorAction*)runManager -> GetUserPrimaryGeneratorAction();
69  // Pointer to the detector material
70  detectorMat = pDet -> GetDetectorLogicalVolume() -> GetMaterial();
71  density = detectorMat -> GetDensity();
72  // Instances for Total LET
73  totalLetD = new G4double[nVoxels];
75 
76 }
77 
79 {
80  Clear();
81  delete [] totalLetD;
82  delete [] DtotalLetD;
83 }
84 
85 // Fill energy spectrum for every voxel (local energy spectrum)
87 {
88  for(G4int v=0; v < nVoxels; v++) totalLetD[v] = DtotalLetD[v] = 0.;
89  Clear();
90 }
95 {
96  for (size_t i=0; i < ionLetStore.size(); i++)
97  {
98  delete [] ionLetStore[i].letDN;
99  delete [] ionLetStore[i].letDD;
100  }
101  ionLetStore.clear();
102 }
104  G4ParticleDefinition* particleDef,
105  /*G4double kinEnergy,*/
106  G4double DE,
107  G4double DX,
108  G4int i, G4int j, G4int k)
109 {
110  if (DE <= 0. || DX <=0.) return;
111  if (!doCalculation) return;
112  G4int Z = particleDef -> GetAtomicNumber();
113 
114 
115  G4int PDGencoding = particleDef -> GetPDGEncoding();
116  PDGencoding -= PDGencoding%10;
117 
118  G4int voxel = matrix -> Index(i,j,k);
119  // Total LET calculation...
120  totalLetD[voxel] += DE*(DE/DX);
121  DtotalLetD[voxel] += DE;
122  // Single ion LET
123  if (Z>=1)
124  {
125  // Search for already allocated data...
126  size_t l;
127  for (l=0; l < ionLetStore.size(); l++)
128  {
129  if (ionLetStore[l].PDGencoding == PDGencoding)
130  if ( ((trackID ==1) && (ionLetStore[l].isPrimary)) || ((trackID !=1) && (!ionLetStore[l].isPrimary)))
131  break;
132  }
133 
134  if (l == ionLetStore.size()) // Just another type of ion/particle for our store...
135  {
136 
137  G4int A = particleDef -> GetAtomicMass();
138 
139  G4String fullName = particleDef -> GetParticleName();
140  G4String name = fullName.substr (0, fullName.find("[") ); // cut excitation energy [x.y]
141 
142  ionLet ion =
143  {
144  (trackID == 1) ? true:false, // is it the primary particle?
145  PDGencoding,
146  fullName,
147  name,
148  Z,
149  A,
150  new G4double[nVoxels], // Let Dose Numerator
151  new G4double[nVoxels] // Let Dose Denominator
152  };
153 
154  // Initialize let
155  for(G4int v=0; v < nVoxels; v++) ion.letDN[v] = ion.letDD[v] = 0.;
156  ionLetStore.push_back(ion);
157  //G4cout << "Allocated LET data for " << ion.name << G4endl;
158 
159  }
160  ionLetStore[l].letDN[voxel] += DE*(DE/DX);
161  ionLetStore[l].letDD[voxel] += DE;
162  }
163 }
164 
165 
166 
167 
169 {
170  for(G4int v=0; v < nVoxels; v++) if (DtotalLetD[v]>0.) totalLetD[v] = totalLetD[v]/DtotalLetD[v];
171  // Sort ions by A and then by Z ...
172  std::sort(ionLetStore.begin(), ionLetStore.end());
173  // Compute Let Track and Let Dose for any single ion
174 
175  for(G4int v=0; v < nVoxels; v++)
176  {
177  for (size_t ion=0; ion < ionLetStore.size(); ion++)
178  {
179  if (ionLetStore[ion].letDD[v] >0.) ionLetStore[ion].letDN[v] = ionLetStore[ion].letDN[v] / ionLetStore[ion].letDD[v];
180 
181  }// end loop over ions
182  }
183 
184 }// end loop over voxels
185 
187 {
188 #define width 15L
189  if(ionLetStore.size())
190  {
191  ofs.open(filename, std::ios::out);
192  if (ofs.is_open())
193  {
194 
195  // Write the voxels index and the list of particles/ions
196  ofs << std::setprecision(6) << std::left <<
197  "i\tj\tk\t";
198  ofs << std::setw(width) << "LDT";
199  for (size_t l=0; l < ionLetStore.size(); l++)
200  {
201  G4String a = (ionLetStore[l].isPrimary) ? "_1":"";
202  ofs << std::setw(width) << ionLetStore[l].name + a ;
203  }
204  ofs << G4endl;
205 
206  // Write data
207  for(G4int i = 0; i < numberOfVoxelAlongX; i++)
208  for(G4int j = 0; j < numberOfVoxelAlongY; j++)
209  for(G4int k = 0; k < numberOfVoxelAlongZ; k++)
210  {
211  G4int v = matrix -> Index(i, j, k);
212  // row write
213  for (size_t l=0; l < ionLetStore.size(); l++)
214  {
215  // Write only not identically null data lines
216  if(ionLetStore[l].letDN)
217  {
218  ofs << G4endl;
219  ofs << i << '\t' << j << '\t' << k << '\t';
220 
221  ofs << std::setw(width) << totalLetD[v]/(keV/um);
222  for (size_t ll=0; ll < ionLetStore.size(); ll++)
223  {
224  ofs << std::setw(width) << ionLetStore[ll].letDN[v]/(keV/um) ;
225  }
226  break;
227  }
228  }
229  }
230  ofs.close();
231  G4cout << "Let is being written to " << filename << G4endl;
232  }
233 
234  }
235 }
236 
238 {
239 #ifdef G4ANALYSIS_USE_ROOT
240 
242 
243  for(G4int i = 0; i < numberOfVoxelAlongX; i++)
244  for(G4int j = 0; j < numberOfVoxelAlongY; j++)
245  for(G4int k = 0; k < numberOfVoxelAlongZ; k++)
246  {
247  G4int v = matrix -> Index(i, j, k);
248  for (size_t ion=0; ion < ionLetStore.size(); ion++)
249  {
250 
251  analysis -> FillLetFragmentTuple( i, j, k, ionLetStore[ion].A, ionLetStore[ion].Z, ionLetStore[ion].letDN[v]);
252 
253 
254  }
255  }
256 
257 #endif
258 }
259 
static G4bool doCalculation
void FillEnergySpectrum(G4int trackID, G4ParticleDefinition *particleDef, G4double DE, G4double DX, G4int i, G4int j, G4int k)
G4double * letDD
static HadrontherapyAnalysisManager * GetInstance()
Get the pointer to the analysis manager.
HadrontherapyMatrix * matrix
static HadrontherapyLet * GetInstance()
HadrontherapyPrimaryGeneratorAction * pPGA
G4String name
Definition: TRTMaterials.hh:40
#define width
G4double a
Definition: TRTMaterials.hh:39
int G4int
Definition: G4Types.hh:78
static HadrontherapyMatrix * GetInstance()
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
HadrontherapyLet(HadrontherapyDetectorConstruction *)
G4Material * detectorMat
std::vector< ionLet > ionLetStore
void Clear()
Clear all stored data.
static const G4double A[nN]
std::ofstream ofs
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static G4RunManager * GetRunManager()
Definition: G4RunManager.cc:79
A class for connecting the simulation to an analysis package.
#define G4endl
Definition: G4ios.hh:61
static HadrontherapyLet * instance
G4double * letDN
static const double keV
Definition: G4SIunits.hh:195
double G4double
Definition: G4Types.hh:76