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