| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // Hadrontherapy advanced example for Geant4
27 // See more at: https://twiki.cern.ch/twiki/bi
28
29 #include "HadrontherapyDetectorConstruction.hh
30 #include "HadrontherapyLet.hh"
31
32 #include "HadrontherapyMatrix.hh"
33 #include "HadrontherapyInteractionParameters.h
34 #include "HadrontherapyPrimaryGeneratorAction.
35 #include "HadrontherapyMatrix.hh"
36 #include "G4AnalysisManager.hh"
37 #include "G4RunManager.hh"
38 #include "G4SystemOfUnits.hh"
39 #include <cmath>
40
41 HadrontherapyLet* HadrontherapyLet::instance =
42 G4bool HadrontherapyLet::doCalculation = false
43
44 HadrontherapyLet* HadrontherapyLet::GetInstanc
45 {
46 if (instance) delete instance;
47 instance = new HadrontherapyLet(pDet);
48 return instance;
49 }
50
51 HadrontherapyLet* HadrontherapyLet::GetInstanc
52 {
53 return instance;
54 }
55
56 HadrontherapyLet::HadrontherapyLet(Hadronthera
57 :filename("Let.out"),matrix(0) // Default outp
58 {
59
60 matrix = HadrontherapyMatrix::GetInstance(
61
62 if (!matrix)
63 G4Exception("HadrontherapyLet::Hadront
64 "Hadrontherapy0005", Fatal
65 "HadrontherapyMatrix not f
66
67 nVoxels = matrix -> GetNvoxel();
68
69 numberOfVoxelAlongX = matrix -> GetNumberO
70 numberOfVoxelAlongY = matrix -> GetNumberO
71 numberOfVoxelAlongZ = matrix -> GetNumberO
72
73 G4RunManager *runManager = G4RunManager::G
74 pPGA = (HadrontherapyPrimaryGeneratorActio
75 // Pointer to the detector material
76 detectorMat = pDet -> GetDetectorLogicalVo
77 density = detectorMat -> GetDensity();
78 // Instances for Total LET
79 totalLetD = new G4double[nVoxels];
80 DtotalLetD = new G4double[nVoxels];
81 totalLetT = new G4double[nVoxels];
82 DtotalLetT = new G4double[nVoxels];
83
84 }
85
86 HadrontherapyLet::~HadrontherapyLet()
87 {
88 Clear();
89 delete [] totalLetD;
90 delete [] DtotalLetD;
91 delete [] totalLetT;
92 delete [] DtotalLetT;
93 }
94
95 // Fill energy spectrum for every voxel (local
96 void HadrontherapyLet::Initialize()
97 {
98 for(G4int v=0; v < nVoxels; v++) totalLetD
99 Clear();
100 }
101 /**
102 * Clear all stored data
103 */
104 void HadrontherapyLet::Clear()
105 {
106 for (size_t i=0; i < ionLetStore.size(); i
107 {
108 delete [] ionLetStore[i].letDN; // Let
109 delete [] ionLetStore[i].letDD; // Let
110 delete [] ionLetStore[i].letTN; // Let
111 delete [] ionLetStore[i].letTD; // Let
112 }
113 ionLetStore.clear();
114 }
115 void HadrontherapyLet::FillEnergySpectrum(G4i
116 G4P
117 G4d
118 con
119 G4d
120 G4d
121 G4d
122 G4i
123 {
124 if (DE <= 0. || DX <=0. ) return; // calcu
125 if (!doCalculation) return;
126
127 // atomic number
128 G4int Z = particleDef -> GetAtomicNumber()
129 if (Z<1) return; // calculate only protons
130
131 G4int PDGencoding = particleDef -> GetPDGE
132 PDGencoding -= PDGencoding%10; // simple P
133
134 G4int voxel = matrix -> Index(i,j,k);
135
136 // ICRU stopping power calculation
137 G4EmCalculator emCal;
138 // use the mean kinetic energy of ions in
139 G4double Lsn = emCal.ComputeElectronicDEDX
140
141
142 // Total LET calculation...
143 totalLetD[voxel] += (DE + DEEletrons) * L
144 DtotalLetD[voxel] += DE + DEEletrons;
145 totalLetT[voxel] += DX * Lsn;
146 DtotalLetT[voxel] += DX;
147
148 // store primary ions and secondary ions
149 size_t l;
150 for (l=0; l < ionLetStore.size(); l++)
151 {
152 // judge species of the current partic
153 if (ionLetStore[l].PDGencoding == PDGe
154 if ( ((trackID ==1) && (ionLetStor
155 break;
156 }
157
158 if (l == ionLetStore.size()) // Just anoth
159 {
160 // mass number
161 G4int A = particleDef -> GetAtomicMass
162
163 // particle name
164 G4String fullName = particleDef -> Get
165 G4String name = fullName.substr (0, fu
166
167 ionLet ion =
168 {
169 (trackID == 1) ? true:false, // is
170 PDGencoding,
171 fullName,
172 name,
173 Z,
174 A,
175 new G4double[nVoxels], // Let Dose
176 new G4double[nVoxels], // Let Dos
177 new G4double[nVoxels], // Let Trac
178 new G4double[nVoxels], // Let Tra
179 };
180
181 // Initialize let and other parameters
182 for(G4int v=0; v < nVoxels; v++)
183 {
184 ion.letDN[v] = ion.letDD[v] = ion.
185 }
186
187
188 ionLetStore.push_back(ion);
189 }
190
191 // calculate ions let and store them
192 ionLetStore[l].letDN[voxel] += (DE + DEEle
193 ionLetStore[l].letDD[voxel] += DE + DEElet
194 ionLetStore[l].letTN[voxel] += DX* Lsn;
195 ionLetStore[l].letTD[voxel] += DX;
196
197 }
198
199
200
201
202 void HadrontherapyLet::LetOutput()
203 {
204 for(G4int v=0; v < nVoxels; v++)
205 {
206 // compute total let
207 if (DtotalLetD[v]>0.) totalLetD[v] = t
208 if (DtotalLetT[v]>0.) totalLetT[v] = t
209 }
210
211 // Sort ions by A and then by Z ...
212 std::sort(ionLetStore.begin(), ionLetStore
213
214
215 // Compute Let Track and Let Dose for ions
216
217 for(G4int v=0; v < nVoxels; v++)
218 {
219
220 for (size_t ion=0; ion < ionLetStore.s
221 {
222 // compute ions let
223 if (ionLetStore[ion].letDD[v] >0.)
224 if (ionLetStore[ion].letTD[v] >0.)
225 } // end loop over ions
226 }
227 } // end loop over voxels
228
229
230
231 void HadrontherapyLet::StoreLetAscii()
232 {
233 #define width 15L
234
235 if(ionLetStore.size())
236 {
237 ofs.open(filename, std::ios::out);
238 if (ofs.is_open())
239 {
240
241 // Write the voxels index and tota
242 ofs << "i" << '\t' << "j" << '\t'
243
244 ofs << '\t' << "LDT";
245 ofs << '\t' << "LTT";
246
247 for (size_t l=0; l < ionLetStore.s
248 {
249 G4String a = (ionLetStore[l].i
250 ofs << '\t' << ionLetStore[l].
251 G4String b = (ionLetStore[l].i
252 ofs << '\t' << ionLetStore[l].
253 }
254
255
256 // Write data
257
258 G4AnalysisManager* LetFragmentTup
259
260 LetFragmentTuple->SetVerboseLevel(
261 LetFragmentTuple->SetFirstHistoId(
262 LetFragmentTuple->SetFirstNtupleId
263 LetFragmentTuple ->OpenFile("Let.c
264
265
266 LetFragmentTuple ->CreateNtuple("c
267
268
269 LetFragmentTuple ->CreateNtupleICo
270 LetFragmentTuple ->CreateNtupleICo
271 LetFragmentTuple ->CreateNtupleICo
272 LetFragmentTuple ->CreateNtupleDCo
273 LetFragmentTuple ->CreateNtupleDCo
274 LetFragmentTuple ->CreateNtupleICo
275 LetFragmentTuple ->CreateNtupleICo
276 LetFragmentTuple ->CreateNtupleDCo
277 LetFragmentTuple ->CreateNtupleDCo
278 LetFragmentTuple ->FinishNtuple();
279
280
281 for(G4int i = 0; i < numberOfVoxel
282 for(G4int j = 0; j < numberOfV
283 for(G4int k = 0; k < numbe
284 {
285 LetFragmentTuple->Fill
286 LetFragmentTuple->Fill
287 LetFragmentTuple->Fill
288
289 G4int v = matrix -> In
290
291 // write total Lets an
292 ofs << G4endl;
293 ofs << i << '\t' << j
294 ofs << '\t' << totalLe
295 ofs << '\t' << totalLe
296
297
298 // write ions Lets
299 for (size_t l=0; l < i
300 {
301
302 // Write only not
303 if(ionLetStore[l].
304 {
305 // write ions
306 ofs << '\t' <<
307 ofs << '\t' <<
308 }
309 }
310
311 LetFragmentTuple->Fill
312 LetFragmentTuple->Fill
313
314
315 for (size_t ll=0; ll <
316 {
317
318
319 LetFragmentTuple->
320 LetFragmentTuple->
321
322
323 LetFragmentTuple->
324 LetFragmentTuple->
325 LetFragmentTuple->
326 }
327 }
328 ofs.close();
329
330 LetFragmentTuple->Write();
331 LetFragmentTuple->CloseFile();
332 }
333
334 }
335
336 }
337
338