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