| 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 // <<
26 // Hadrontherapy advanced example for Geant4 <<
27 // See more at: https://twiki.cern.ch/twiki/bi <<
28 <<
29 #include "HadrontherapyDetectorConstruction.hh 25 #include "HadrontherapyDetectorConstruction.hh"
30 #include "HadrontherapyLet.hh" 26 #include "HadrontherapyLet.hh"
31 << 27 #include "HadrontherapyAnalysisManager.hh"
32 #include "HadrontherapyMatrix.hh" <<
33 #include "HadrontherapyInteractionParameters.h 28 #include "HadrontherapyInteractionParameters.hh"
34 #include "HadrontherapyPrimaryGeneratorAction. 29 #include "HadrontherapyPrimaryGeneratorAction.hh"
35 #include "HadrontherapyMatrix.hh" 30 #include "HadrontherapyMatrix.hh"
36 #include "G4AnalysisManager.hh" <<
37 #include "G4RunManager.hh" 31 #include "G4RunManager.hh"
38 #include "G4SystemOfUnits.hh" 32 #include "G4SystemOfUnits.hh"
39 #include <cmath> 33 #include <cmath>
40 34
41 HadrontherapyLet* HadrontherapyLet::instance = 35 HadrontherapyLet* HadrontherapyLet::instance = NULL;
42 G4bool HadrontherapyLet::doCalculation = false 36 G4bool HadrontherapyLet::doCalculation = false;
43 37
44 HadrontherapyLet* HadrontherapyLet::GetInstanc 38 HadrontherapyLet* HadrontherapyLet::GetInstance(HadrontherapyDetectorConstruction *pDet)
45 { 39 {
46 if (instance) delete instance; << 40 if (instance) delete instance;
47 instance = new HadrontherapyLet(pDet); << 41 instance = new HadrontherapyLet(pDet);
48 return instance; << 42 return instance;
49 } 43 }
50 44
51 HadrontherapyLet* HadrontherapyLet::GetInstanc 45 HadrontherapyLet* HadrontherapyLet::GetInstance()
52 { 46 {
53 return instance; << 47 return instance;
54 } 48 }
55 49
56 HadrontherapyLet::HadrontherapyLet(Hadronthera 50 HadrontherapyLet::HadrontherapyLet(HadrontherapyDetectorConstruction* pDet)
57 :filename("Let.out"),matrix(0) // Default outp << 51 :filename("Let.out"),matrix(0) // Default output filename
58 { 52 {
59 53
60 matrix = HadrontherapyMatrix::GetInstance( 54 matrix = HadrontherapyMatrix::GetInstance();
61 55
62 if (!matrix) << 56 if (!matrix)
63 G4Exception("HadrontherapyLet::Hadront << 57 G4Exception("HadrontherapyLet::HadrontherapyLet",
64 "Hadrontherapy0005", Fatal << 58 "Hadrontherapy0005", FatalException,
65 "HadrontherapyMatrix not f << 59 "HadrontherapyMatrix not found. Firstly create an instance of it.");
66 60
67 nVoxels = matrix -> GetNvoxel(); 61 nVoxels = matrix -> GetNvoxel();
68 62
69 numberOfVoxelAlongX = matrix -> GetNumberO 63 numberOfVoxelAlongX = matrix -> GetNumberOfVoxelAlongX();
70 numberOfVoxelAlongY = matrix -> GetNumberO 64 numberOfVoxelAlongY = matrix -> GetNumberOfVoxelAlongY();
71 numberOfVoxelAlongZ = matrix -> GetNumberO 65 numberOfVoxelAlongZ = matrix -> GetNumberOfVoxelAlongZ();
72 66
73 G4RunManager *runManager = G4RunManager::G 67 G4RunManager *runManager = G4RunManager::GetRunManager();
74 pPGA = (HadrontherapyPrimaryGeneratorActio 68 pPGA = (HadrontherapyPrimaryGeneratorAction*)runManager -> GetUserPrimaryGeneratorAction();
75 // Pointer to the detector material 69 // Pointer to the detector material
76 detectorMat = pDet -> GetDetectorLogicalVo 70 detectorMat = pDet -> GetDetectorLogicalVolume() -> GetMaterial();
77 density = detectorMat -> GetDensity(); 71 density = detectorMat -> GetDensity();
78 // Instances for Total LET 72 // Instances for Total LET
79 totalLetD = new G4double[nVoxels]; 73 totalLetD = new G4double[nVoxels];
80 DtotalLetD = new G4double[nVoxels]; 74 DtotalLetD = new G4double[nVoxels];
81 totalLetT = new G4double[nVoxels]; <<
82 DtotalLetT = new G4double[nVoxels]; <<
83 75
84 } 76 }
85 77
86 HadrontherapyLet::~HadrontherapyLet() 78 HadrontherapyLet::~HadrontherapyLet()
87 { 79 {
88 Clear(); << 80 Clear();
89 delete [] totalLetD; << 81 delete [] totalLetD;
90 delete [] DtotalLetD; << 82 delete [] DtotalLetD;
91 delete [] totalLetT; <<
92 delete [] DtotalLetT; <<
93 } 83 }
94 84
95 // Fill energy spectrum for every voxel (local 85 // Fill energy spectrum for every voxel (local energy spectrum)
96 void HadrontherapyLet::Initialize() 86 void HadrontherapyLet::Initialize()
97 { 87 {
98 for(G4int v=0; v < nVoxels; v++) totalLetD << 88 for(G4int v=0; v < nVoxels; v++) totalLetD[v] = DtotalLetD[v] = 0.;
99 Clear(); << 89 Clear();
100 } 90 }
101 /** 91 /**
102 * Clear all stored data 92 * Clear all stored data
103 */ 93 */
104 void HadrontherapyLet::Clear() 94 void HadrontherapyLet::Clear()
105 { 95 {
106 for (size_t i=0; i < ionLetStore.size(); i << 96 for (size_t i=0; i < ionLetStore.size(); i++)
107 { << 97 {
108 delete [] ionLetStore[i].letDN; // Let << 98 delete [] ionLetStore[i].letDN;
109 delete [] ionLetStore[i].letDD; // Let << 99 delete [] ionLetStore[i].letDD;
110 delete [] ionLetStore[i].letTN; // Let << 100 }
111 delete [] ionLetStore[i].letTD; // Let << 101 ionLetStore.clear();
112 } <<
113 ionLetStore.clear(); <<
114 } 102 }
115 void HadrontherapyLet::FillEnergySpectrum(G4i 103 void HadrontherapyLet::FillEnergySpectrum(G4int trackID,
116 G4P 104 G4ParticleDefinition* particleDef,
117 G4d << 105 /*G4double kinEnergy,*/
118 con <<
119 G4d 106 G4double DE,
120 G4d <<
121 G4d 107 G4double DX,
122 G4i 108 G4int i, G4int j, G4int k)
123 { 109 {
124 if (DE <= 0. || DX <=0. ) return; // calcu << 110 if (DE <= 0. || DX <=0.) return;
125 if (!doCalculation) return; << 111 if (!doCalculation) return;
126 << 112 G4int Z = particleDef -> GetAtomicNumber();
127 // atomic number << 113
128 G4int Z = particleDef -> GetAtomicNumber() << 114
129 if (Z<1) return; // calculate only protons << 115 G4int PDGencoding = particleDef -> GetPDGEncoding();
130 << 116 PDGencoding -= PDGencoding%10;
131 G4int PDGencoding = particleDef -> GetPDGE << 117
132 PDGencoding -= PDGencoding%10; // simple P << 118 G4int voxel = matrix -> Index(i,j,k);
133 << 119 // Total LET calculation...
134 G4int voxel = matrix -> Index(i,j,k); << 120 totalLetD[voxel] += DE*(DE/DX);
135 << 121 DtotalLetD[voxel] += DE;
136 // ICRU stopping power calculation << 122 // Single ion LET
137 G4EmCalculator emCal; << 123 if (Z>=1)
138 // use the mean kinetic energy of ions in << 124 {
139 G4double Lsn = emCal.ComputeElectronicDEDX << 125 // Search for already allocated data...
140 << 126 size_t l;
141 << 127 for (l=0; l < ionLetStore.size(); l++)
142 // Total LET calculation... << 128 {
143 totalLetD[voxel] += (DE + DEEletrons) * L << 129 if (ionLetStore[l].PDGencoding == PDGencoding)
144 DtotalLetD[voxel] += DE + DEEletrons; << 130 if ( ((trackID ==1) && (ionLetStore[l].isPrimary)) || ((trackID !=1) && (!ionLetStore[l].isPrimary)))
145 totalLetT[voxel] += DX * Lsn; << 131 break;
146 DtotalLetT[voxel] += DX; << 132 }
147 << 133
148 // store primary ions and secondary ions << 134 if (l == ionLetStore.size()) // Just another type of ion/particle for our store...
149 size_t l; << 135 {
150 for (l=0; l < ionLetStore.size(); l++) << 136
151 { << 137 G4int A = particleDef -> GetAtomicMass();
152 // judge species of the current partic << 138
153 if (ionLetStore[l].PDGencoding == PDGe << 139 G4String fullName = particleDef -> GetParticleName();
154 if ( ((trackID ==1) && (ionLetStor << 140 G4String name = fullName.substr (0, fullName.find("[") ); // cut excitation energy [x.y]
155 break; << 141
156 } << 142 ionLet ion =
157 << 143 {
158 if (l == ionLetStore.size()) // Just anoth << 144 (trackID == 1) ? true:false, // is it the primary particle?
159 { << 145 PDGencoding,
160 // mass number << 146 fullName,
161 G4int A = particleDef -> GetAtomicMass << 147 name,
162 << 148 Z,
163 // particle name << 149 A,
164 G4String fullName = particleDef -> Get << 150 new G4double[nVoxels], // Let Dose Numerator
165 G4String name = fullName.substr (0, fu << 151 new G4double[nVoxels] // Let Dose Denominator
166 << 152 };
167 ionLet ion = << 153
168 { << 154 // Initialize let
169 (trackID == 1) ? true:false, // is << 155 for(G4int v=0; v < nVoxels; v++) ion.letDN[v] = ion.letDD[v] = 0.;
170 PDGencoding, << 156 ionLetStore.push_back(ion);
171 fullName, << 157 //G4cout << "Allocated LET data for " << ion.name << G4endl;
172 name, << 158
173 Z, << 159 }
174 A, << 160 ionLetStore[l].letDN[voxel] += DE*(DE/DX);
175 new G4double[nVoxels], // Let Dose << 161 ionLetStore[l].letDD[voxel] += DE;
176 new G4double[nVoxels], // Let Dos << 162 }
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 } 163 }
198 164
199 165
200 166
201 167
202 void HadrontherapyLet::LetOutput() 168 void HadrontherapyLet::LetOutput()
203 { 169 {
204 for(G4int v=0; v < nVoxels; v++) << 170 for(G4int v=0; v < nVoxels; v++) if (DtotalLetD[v]>0.) totalLetD[v] = totalLetD[v]/DtotalLetD[v];
205 { << 171 // Sort ions by A and then by Z ...
206 // compute total let << 172 std::sort(ionLetStore.begin(), ionLetStore.end());
207 if (DtotalLetD[v]>0.) totalLetD[v] = t << 173 // Compute Let Track and Let Dose for any single ion
208 if (DtotalLetT[v]>0.) totalLetT[v] = t << 174
209 } << 175 for(G4int v=0; v < nVoxels; v++)
210 << 176 {
211 // Sort ions by A and then by Z ... << 177 for (size_t ion=0; ion < ionLetStore.size(); ion++)
212 std::sort(ionLetStore.begin(), ionLetStore << 178 {
213 << 179 if (ionLetStore[ion].letDD[v] >0.) ionLetStore[ion].letDN[v] = ionLetStore[ion].letDN[v] / ionLetStore[ion].letDD[v];
214 << 180
215 // Compute Let Track and Let Dose for ions << 181 }// end loop over ions
>> 182 }
216 183
217 for(G4int v=0; v < nVoxels; v++) << 184 }// end loop over voxels
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 185
231 void HadrontherapyLet::StoreLetAscii() 186 void HadrontherapyLet::StoreLetAscii()
232 { 187 {
233 #define width 15L 188 #define width 15L
234 << 189 if(ionLetStore.size())
235 if(ionLetStore.size()) << 190 {
236 { << 191 ofs.open(filename, std::ios::out);
237 ofs.open(filename, std::ios::out); << 192 if (ofs.is_open())
238 if (ofs.is_open()) << 193 {
239 { << 194
240 << 195 // Write the voxels index and the list of particles/ions
241 // Write the voxels index and tota << 196 ofs << std::setprecision(6) << std::left <<
242 ofs << "i" << '\t' << "j" << '\t' << 197 "i\tj\tk\t";
243 << 198 ofs << std::setw(width) << "LDT";
244 ofs << '\t' << "LDT"; << 199 for (size_t l=0; l < ionLetStore.size(); l++)
245 ofs << '\t' << "LTT"; << 200 {
246 << 201 G4String a = (ionLetStore[l].isPrimary) ? "_1":"";
247 for (size_t l=0; l < ionLetStore.s << 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
>> 237 void HadrontherapyLet::StoreLetRoot()
>> 238 {
>> 239 #ifdef G4ANALYSIS_USE_ROOT
>> 240
>> 241 HadrontherapyAnalysisManager* analysis = HadrontherapyAnalysisManager::GetInstance();
>> 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++)
248 { 246 {
249 G4String a = (ionLetStore[l].i << 247 G4int v = matrix -> Index(i, j, k);
250 ofs << '\t' << ionLetStore[l]. << 248 for (size_t ion=0; ion < ionLetStore.size(); ion++)
251 G4String b = (ionLetStore[l].i << 249 {
252 ofs << '\t' << ionLetStore[l]. << 250
>> 251 analysis -> FillLetFragmentTuple( i, j, k, ionLetStore[ion].A, ionLetStore[ion].Z, ionLetStore[ion].letDN[v]);
>> 252
>> 253
>> 254 }
253 } 255 }
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 256
>> 257 #endif
336 } 258 }
337 259
338 260