| 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 // This is the *BASIC* version of Hadrontherapy, a Geant4-based application
27 // See more at: https://twiki.cern.ch/twiki/bi << 27 // See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
>> 28 //
>> 29 // Visit the Hadrontherapy web site (http://www.lns.infn.it/link/Hadrontherapy) to request
>> 30 // the *COMPLETE* version of this program, together with its documentation;
>> 31 // Hadrontherapy (both basic and full version) are supported by the Italian INFN
>> 32 // Institute in the framework of the MC-INFN Group
>> 33 //
>> 34
>> 35 #include "HadrontherapyMatrix.hh"
>> 36 #include "HadrontherapyAnalysisManager.hh"
>> 37 #include "G4RunManager.hh"
>> 38 #include "HadrontherapyPrimaryGeneratorAction.hh"
>> 39 #include "G4ParticleGun.hh"
28 40
29 #include <fstream> 41 #include <fstream>
30 #include <iostream> 42 #include <iostream>
31 #include <sstream> 43 #include <sstream>
32 #include <iomanip> 44 #include <iomanip>
33 <<
34 #include "HadrontherapyMatrix.hh" <<
35 #include "HadrontherapyPrimaryGeneratorAction. <<
36 #include "globals.hh" 45 #include "globals.hh"
37 #include "G4SystemOfUnits.hh" <<
38 #include "G4RunManager.hh" <<
39 #include "G4ParticleGun.hh" <<
40 #include "HadrontherapySteppingAction.hh" <<
41 #include "HadrontherapyAnalysisFileMessenger.h <<
42 #include "G4SystemOfUnits.hh" <<
43 #include <time.h> <<
44 <<
45 HadrontherapyAnalysis* HadrontherapyAnalysis:: <<
46 ////////////////////////////////////////////// <<
47 <<
48 HadrontherapyAnalysis::HadrontherapyAnalysis() <<
49 { <<
50 fMess = new HadrontherapyAnalysisFileMesse <<
51 } <<
52 <<
53 ////////////////////////////////////////////// <<
54 HadrontherapyAnalysis::~HadrontherapyAnalysis( <<
55 { <<
56 delete fMess; <<
57 } <<
58 <<
59 ////////////////////////////////////////////// <<
60 HadrontherapyAnalysis* HadrontherapyAnalysis:: <<
61 <<
62 if (instance == 0) instance = new Hadronth <<
63 return instance; <<
64 } <<
65 46
>> 47 // Units definition: CLHEP/Units/SystemOfUnits.h
>> 48 //
66 HadrontherapyMatrix* HadrontherapyMatrix::inst 49 HadrontherapyMatrix* HadrontherapyMatrix::instance = NULL;
67 G4bool HadrontherapyMatrix::secondary = false; 50 G4bool HadrontherapyMatrix::secondary = false;
68 51
69 <<
70 // Only return a pointer to matrix 52 // Only return a pointer to matrix
71 HadrontherapyMatrix* HadrontherapyMatrix::GetI << 53 HadrontherapyMatrix* HadrontherapyMatrix::GetInstance()
72 { 54 {
73 return instance; << 55 return instance;
74 } 56 }
75 << 57 // This STATIC method delete (!) the old matrix and rewrite a new object returning a pointer to it
76 ////////////////////////////////////////////// << 58 // TODO A check on the parameters is required!
77 // This STATIC method delete (!) the old matri << 59 HadrontherapyMatrix* HadrontherapyMatrix::GetInstance(G4int voxelX, G4int voxelY, G4int voxelZ, G4double mass)
78 // TODO A check on the parameters is required! << 60 {
79 HadrontherapyMatrix* HadrontherapyMatrix::GetI << 61 if (instance) delete instance;
80 { << 62 instance = new HadrontherapyMatrix(voxelX, voxelY, voxelZ, mass);
81 if (instance) delete instance; << 63 instance -> Initialize();
82 instance = new HadrontherapyMatrix(voxelX, << 64 return instance;
83 instance -> Initialize(); <<
84 return instance; <<
85 } 65 }
86 <<
87 ////////////////////////////////////////////// <<
88 HadrontherapyMatrix::HadrontherapyMatrix(G4int 66 HadrontherapyMatrix::HadrontherapyMatrix(G4int voxelX, G4int voxelY, G4int voxelZ, G4double mass):
89 stdFile("Dose.out"), 67 stdFile("Dose.out"),
90 doseUnit(gray) << 68 doseUnit(MeV/g)
91 { << 69 {
92 // Number of the voxels of the phantom << 70 // Number of the voxels of the phantom
93 // For Y = Z = 1 the phantom is divided in << 71 // For Y = Z = 1 the phantom is divided in slices (and not in voxels)
94 // orthogonal to the beam axis << 72 // orthogonal to the beam axis
95 numberOfVoxelAlongX = voxelX; << 73 numberOfVoxelAlongX = voxelX;
96 numberOfVoxelAlongY = voxelY; << 74 numberOfVoxelAlongY = voxelY;
97 numberOfVoxelAlongZ = voxelZ; << 75 numberOfVoxelAlongZ = voxelZ;
98 massOfVoxel = mass; << 76 massOfVoxel = mass;
99 << 77 // Create the dose matrix
100 << 78 matrix = new G4double[numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ];
101 // Create the dose matrix << 79 if (matrix)
102 matrix = new G4double[numberOfVoxelAlongX* << 80 {
103 if (matrix) << 81 G4cout << "HadrontherapyMatrix: Memory space to store physical dose into " <<
104 { << 82 numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ <<
105 G4cout << "HadrontherapyMatrix: Memory << 83 " voxels has been allocated " << G4endl;
106 numberOfVoxelAlongX*numberOf << 84 }
107 " voxels has been allocated << 85 else G4Exception("HadrontherapyMatrix::HadrontherapyMatrix()", "Hadrontherapy0005", FatalException, "Can't allocate memory to store physical dose!");
108 } << 86 // Hit voxel (TrackID) marker
109 << 87 // This array mark the status of voxel, if a hit occur, with the trackID of the particle
110 else G4Exception("HadrontherapyMatrix::Had << 88 // Must be initialized
111 << 89 hitTrack = new G4int[numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ];
112 << 90 ClearHitTrack();
113 // Hit voxel (TrackID) marker <<
114 // This array mark the status of voxel, if <<
115 // Must be initialized <<
116 <<
117 hitTrack = new G4int[numberOfVoxelAlongX*n <<
118 ClearHitTrack(); <<
119 } 91 }
120 92
121 ////////////////////////////////////////////// 93 /////////////////////////////////////////////////////////////////////////////
122 HadrontherapyMatrix::~HadrontherapyMatrix() 94 HadrontherapyMatrix::~HadrontherapyMatrix()
123 { 95 {
124 delete[] matrix; 96 delete[] matrix;
125 delete[] hitTrack; 97 delete[] hitTrack;
>> 98 // free fluences/dose data memory
126 Clear(); 99 Clear();
127 } 100 }
128 101
129 ////////////////////////////////////////////// 102 /////////////////////////////////////////////////////////////////////////////
130 void HadrontherapyMatrix::Clear() 103 void HadrontherapyMatrix::Clear()
131 { 104 {
132 for (size_t i=0; i<ionStore.size(); i++) 105 for (size_t i=0; i<ionStore.size(); i++)
133 { 106 {
134 delete[] ionStore[i].dose; << 107 delete[] ionStore[i].dose;
135 delete[] ionStore[i].fluence; << 108 delete[] ionStore[i].fluence;
136 } 109 }
137 ionStore.clear(); 110 ionStore.clear();
138 } 111 }
139 112
140 ////////////////////////////////////////////// 113 /////////////////////////////////////////////////////////////////////////////
141 // Initialise the elements of the matrix to ze 114 // Initialise the elements of the matrix to zero
142 <<
143 void HadrontherapyMatrix::Initialize() 115 void HadrontherapyMatrix::Initialize()
144 { << 116 {
145 // Clear ions store 117 // Clear ions store
146 Clear(); 118 Clear();
147 // Clear dose 119 // Clear dose
148 for(int i=0;i<numberOfVoxelAlongX*numberOf 120 for(int i=0;i<numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ;i++)
149 { 121 {
150 matrix[i] = 0; << 122 matrix[i] = 0;
151 } 123 }
152 } 124 }
153 << 125 /////////////////////////////////////////////////////////////////////////////
154 ////////////////////////////////////////////// << 126 /////////////////////////////////////////////////////////////////////////////
155 // Print generated nuclides list << 127 // Print generated nuclides list
156 <<
157 ////////////////////////////////////////////// <<
158 void HadrontherapyMatrix::PrintNuclides() 128 void HadrontherapyMatrix::PrintNuclides()
159 { 129 {
160 for (size_t i=0; i<ionStore.size(); i++) 130 for (size_t i=0; i<ionStore.size(); i++)
161 { 131 {
162 G4cout << ionStore[i].name << G4endl; << 132 G4cout << ionStore[i].name << G4endl;
163 } 133 }
164 } 134 }
165 << 135 /////////////////////////////////////////////////////////////////////////////
166 ////////////////////////////////////////////// << 136 // Clear Hit voxel (TrackID) markers
167 // Clear Hit voxel (TrackID) markers <<
168 <<
169 void HadrontherapyMatrix::ClearHitTrack() 137 void HadrontherapyMatrix::ClearHitTrack()
170 { 138 {
171 for(G4int i=0; i<numberOfVoxelAlongX*numbe << 139 for(G4int i=0; i<numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ; i++) hitTrack[i] = 0;
172 } 140 }
173 << 141 // Return Hit status
174 // Return Hit status <<
175 G4int* HadrontherapyMatrix::GetHitTrack(G4int 142 G4int* HadrontherapyMatrix::GetHitTrack(G4int i, G4int j, G4int k)
176 { 143 {
177 return &(hitTrack[Index(i,j,k)]); << 144 return &(hitTrack[Index(i,j,k)]);
178 } 145 }
179 <<
180 ////////////////////////////////////////////// 146 /////////////////////////////////////////////////////////////////////////////
181 // Dose methods... 147 // Dose methods...
182 // Fill DOSE/fluence matrix for secondary part << 148 // Fill DOSE/fluence matrix for secondary particles:
183 // If fluence parameter is true (default value << 149 // If fluence parameter is true (default value is FALSE) then fluence at voxel (i, j, k) is increased.
184 // The energyDeposit parameter fill the dose m << 150 // The energyDeposit parameter fill the dose matrix for voxel (i,j,k)
185 ////////////////////////////////////////////// 151 /////////////////////////////////////////////////////////////////////////////
>> 152
186 G4bool HadrontherapyMatrix::Fill(G4int trackID 153 G4bool HadrontherapyMatrix::Fill(G4int trackID,
187 G4ParticleDef << 154 G4ParticleDefinition* particleDef,
188 G4int i, G4in << 155 G4int i, G4int j, G4int k,
189 G4double ener << 156 G4double energyDeposit,
190 G4bool fluenc << 157 G4bool fluence)
191 { 158 {
192 <<
193 if ( (energyDeposit <=0. && !fluence) || ! 159 if ( (energyDeposit <=0. && !fluence) || !secondary) return false;
194 << 160 // Get Particle Data Group particle ID
195 // Get Particle Data Group particle ID <<
196 G4int PDGencoding = particleDef -> GetPDGE 161 G4int PDGencoding = particleDef -> GetPDGEncoding();
197 PDGencoding -= PDGencoding%10; 162 PDGencoding -= PDGencoding%10;
198 << 163
199 // Search for already allocated data... 164 // Search for already allocated data...
200 for (size_t l=0; l < ionStore.size(); l++) 165 for (size_t l=0; l < ionStore.size(); l++)
201 { 166 {
202 if (ionStore[l].PDGencoding == PDGenco << 167 if (ionStore[l].PDGencoding == PDGencoding )
203 { // Is it a primary or a secondary << 168 { // Is it a primary or a secondary particle?
204 << 169 if ( (trackID ==1 && ionStore[l].isPrimary) || (trackID !=1 && !ionStore[l].isPrimary))
205 if ( (trackID ==1 && ionStore[l].i << 170 {
206 { << 171 if (energyDeposit > 0.) ionStore[l].dose[Index(i, j, k)] += energyDeposit;
207 if (energyDeposit > 0.) << 172
208 << 173 // Fill a matrix per each ion with the fluence
209 ionStore[l].dose[Index(i, << 174 if (fluence) ionStore[l].fluence[Index(i, j, k)]++;
210 << 175 return true;
211 // Fill a matrix per each ion << 176 }
212 << 177 }
213 if (fluence) ionStore[l].fluen <<
214 return true; <<
215 } <<
216 } <<
217 } 178 }
>> 179
218 G4int Z = particleDef-> GetAtomicNumber(); 180 G4int Z = particleDef-> GetAtomicNumber();
219 G4int A = particleDef-> GetAtomicMass(); 181 G4int A = particleDef-> GetAtomicMass();
220 G4String fullName = particleDef -> GetPart <<
221 G4String name = fullName.substr (0, fullNa <<
222 182
223 // Let's put a new particle in our store.. << 183 G4String fullName = particleDef -> GetParticleName();
224 ion newIon = << 184 G4String name = fullName.substr (0, fullName.find("[") ); // cut excitation energy
>> 185 // Let's put a new particle in our store...
>> 186 ion newIon =
225 { 187 {
226 (trackID == 1) ? true:false, << 188 (trackID == 1) ? true:false,
227 PDGencoding, << 189 PDGencoding,
228 name, << 190 name,
229 name.length(), << 191 name.length(),
230 Z, << 192 Z,
231 A, << 193 A,
232 new G4double[numberOfVoxelAlongX * num << 194 new G4double[numberOfVoxelAlongX * numberOfVoxelAlongY * numberOfVoxelAlongZ],
233 new unsigned int[numberOfVoxelAlongX * << 195 new unsigned int[numberOfVoxelAlongX * numberOfVoxelAlongY * numberOfVoxelAlongZ]
234 }; << 196 };
235 << 197 // Initialize data
236 <<
237 // Initialize data <<
238 if (newIon.dose && newIon.fluence) 198 if (newIon.dose && newIon.fluence)
239 { 199 {
240 for(G4int q=0; q<numberOfVoxelAlongX*n << 200 for(G4int m=0; m<numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ; m++)
241 { << 201 {
242 newIon.dose[q] = 0.; << 202 newIon.dose[m] = 0.;
243 newIon.fluence[q] = 0; << 203 newIon.fluence[m] = 0;
244 } << 204 }
245 << 205 if (energyDeposit > 0.) newIon.dose[Index(i, j, k)] += energyDeposit;
246 if (energyDeposit > 0.) newIon.dose[In << 206 if (fluence) newIon.fluence[Index(i, j, k)]++;
247 if (fluence) newIon.fluence[Index(i, j << 207
248 << 208 ionStore.push_back(newIon);
249 ionStore.push_back(newIon); << 209
250 return true; << 210 // TODO Put some verbosity check
>> 211 /*
>> 212 G4cout << "Memory space to store the DOSE/FLUENCE into " <<
>> 213 numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ <<
>> 214 " voxels has been allocated for the nuclide " << newIon.name <<
>> 215 " (Z = " << Z << ", A = " << A << ")" << G4endl ;
>> 216 */
>> 217 return true;
251 } 218 }
252 <<
253 else // XXX Out of memory! XXX 219 else // XXX Out of memory! XXX
254 <<
255 { 220 {
256 return false; << 221 return false;
257 } 222 }
>> 223
258 } 224 }
259 << 225
260 ////////////////////////////////////////////// << 226 /////////////////////////////////////////////////////////////////////////////
261 ////////////////////////////////////////////// << 227 /////////////////////////////////////////////////////////////////////////////
262 // Methods to store data to filenames... << 228 // Methods to store data to filenames...
263 ////////////////////////////////////////////// << 229 ////////////////////////////////////////////////////////////////////////////
264 ////////////////////////////////////////////// << 230 ////////////////////////////////////////////////////////////////////////////
265 // << 231 //
266 // General method to store matrix data to file << 232 // General method to store matrix data to filename
267 void HadrontherapyMatrix::StoreMatrix(G4String 233 void HadrontherapyMatrix::StoreMatrix(G4String file, void* data, size_t psize)
268 { 234 {
269 if (data) 235 if (data)
270 { 236 {
271 ofs.open(file, std::ios::out); << 237 ofs.open(file, std::ios::out);
272 if (ofs.is_open()) << 238 if (ofs.is_open())
273 { << 239 {
274 for(G4int i = 0; i < numberOfVoxel << 240 for(G4int i = 0; i < numberOfVoxelAlongX; i++)
275 for(G4int j = 0; j < numberOfV << 241 for(G4int j = 0; j < numberOfVoxelAlongY; j++)
276 for(G4int k = 0; k < numbe << 242 for(G4int k = 0; k < numberOfVoxelAlongZ; k++)
277 { << 243 {
278 G4int n = Index(i, j, << 244 G4int n = Index(i, j, k);
279 << 245 // Check for data type: u_int, G4double, XXX
280 if (psize == sizeof(un << 246 if (psize == sizeof(unsigned int))
281 { << 247 {
282 unsigned int* pdat << 248 unsigned int* pdata = (unsigned int*)data;
283 << 249 if (pdata[n]) ofs << i << '\t' << j << '\t' <<
284 if (pdata[n]) << 250 k << '\t' << pdata[n] << G4endl;
285 << 251 }
286 ofs << i << '\ << 252 else if (psize == sizeof(G4double))
287 } << 253 {
288 << 254 G4double* pdata = (G4double*)data;
289 else if (psize == size << 255 if (pdata[n]) ofs << i << '\t' << j << '\t' <<
290 << 256 k << '\t' << pdata[n] << G4endl;
291 { << 257 }
292 G4double* pdata = << 258 }
293 if (pdata[n]) ofs << 259 ofs.close();
294 } << 260 }
295 } <<
296 ofs.close(); <<
297 } <<
298 } 261 }
299 } 262 }
300 263
301 ////////////////////////////////////////////// << 264 // Store fluence per single ion in multiple files
302 // Store fluence per single ion in multiple fi <<
303 void HadrontherapyMatrix::StoreFluenceData() 265 void HadrontherapyMatrix::StoreFluenceData()
304 { 266 {
305 for (size_t i=0; i < ionStore.size(); i++) 267 for (size_t i=0; i < ionStore.size(); i++){
306 StoreMatrix(ionStore[i].name + "_Fluen << 268 StoreMatrix(ionStore[i].name + "_Fluence.out", ionStore[i].fluence, sizeof(unsigned int));
307 } 269 }
308 } 270 }
309 << 271 // Store dose per single ion in multiple files
310 ////////////////////////////////////////////// <<
311 // Store dose per single ion in multiple files <<
312 void HadrontherapyMatrix::StoreDoseData() 272 void HadrontherapyMatrix::StoreDoseData()
313 { 273 {
314 << 274
315 for (size_t i=0; i < ionStore.size(); i++) 275 for (size_t i=0; i < ionStore.size(); i++){
316 StoreMatrix(ionStore[i].name + "_Dose. << 276 StoreMatrix(ionStore[i].name + "_Dose.out", ionStore[i].dose, sizeof(G4double));
317 } 277 }
318 } 278 }
319 << 279 /////////////////////////////////////////////////////////////////////////
320 ////////////////////////////////////////////// << 280 // Store dose for all ions into a single file and into ntuples.
321 // Store dose into a single file << 281 // Please note that this function is called via messenger commands
322 // or in histograms. Please, note that this fu << 282 // defined in the HadrontherapyAnalysisFileMessenger.cc class file
323 // messenger commands <<
324 // defined in the HadrontherapyAnalysisFileMes <<
325 void HadrontherapyMatrix::StoreDoseFluenceAsci 283 void HadrontherapyMatrix::StoreDoseFluenceAscii(G4String file)
326 { 284 {
327 #define width 15L 285 #define width 15L
328 filename = (file=="") ? stdFile:file; 286 filename = (file=="") ? stdFile:file;
329 <<
330 // Sort like periodic table 287 // Sort like periodic table
331 <<
332 std::sort(ionStore.begin(), ionStore.end() 288 std::sort(ionStore.begin(), ionStore.end());
333 G4cout << "Dose is being written to " << f 289 G4cout << "Dose is being written to " << filename << G4endl;
334 ofs.open(filename, std::ios::out); 290 ofs.open(filename, std::ios::out);
335 <<
336 if (ofs.is_open()) 291 if (ofs.is_open())
337 { 292 {
338 // Write the voxels index and the list << 293 // Write the voxels index and the list of particles/ions
339 //ofs << std::setprecision(6) << std:: << 294 ofs << std::setprecision(6) << std::left <<
340 ofs << "i" << '\t' << "j" << '\t' << " << 295 "i\tj\tk\t";
341 G4cout << "i" << '\t' << "j" << '\t' < << 296 // Total dose
342 << 297 ofs << std::setw(width) << "Dose(MeV/g)";
343 // Total dose << 298 if (secondary)
344 ofs <<'\t' <<"Dose(Gy)"; << 299 {
345 //ofs << std::setw(width) << "Dose(Gy) << 300 for (size_t l=0; l < ionStore.size(); l++)
346 G4cout << '\t' << "Dose(Gy)"; << 301 {
347 << 302 G4String a = (ionStore[l].isPrimary) ? "_1":""; // is it a primary?
348 G4String fluence = "_f"; << 303 ofs << std::setw(width) << ionStore[l].name + a <<
349 if (secondary) << 304 std::setw(width) << ionStore[l].name + a;
350 { << 305 }
351 for (size_t l=0; l < ionStore.size << 306 ofs << G4endl;
352 { << 307
353 G4String a = (ionStore[l].isPr << 308 /*
354 << 309 * PDGencondig
355 // ofs << std::setw(width) << i << 310 */
356 // std::setw(width) << i << 311 /*
357 << 312 ofs << std::setprecision(6) << std::left <<
358 ofs << '\t' << ionStore[l].nam << 313 "0\t0\t0\t";
359 '\t' << ionStore[l]. << 314
360 << 315 // Total dose
361 G4cout << '\t' << ionStore[l]. << 316 ofs << std::setw(width) << '0';
362 '\t' << ionStore[l]. << 317 for (size_t l=0; l < ionStore.size(); l++)
363 << 318 {
364 } << 319 ofs << std::setw(width) << ionStore[l].PDGencoding <<
365 //ofs << G4endl; << 320 std::setw(width) << ionStore[l].PDGencoding;
366 } << 321 }
367 << 322 ofs << G4endl;
368 // Write data << 323 */
369 for(G4int i = 0; i < numberOfVoxelAlon << 324 }
370 for(G4int j = 0; j < numberOfVoxel << 325 // Write data
371 for(G4int k = 0; k < numberOfV << 326 for(G4int i = 0; i < numberOfVoxelAlongX; i++)
372 { << 327 for(G4int j = 0; j < numberOfVoxelAlongY; j++)
373 G4int n = Index(i, j, k); << 328 for(G4int k = 0; k < numberOfVoxelAlongZ; k++)
374 << 329 {
375 if (matrix[n]) << 330 G4int n = Index(i, j, k);
376 { << 331 // Write only not identically null data lines
377 ofs << G4endl; << 332 if (matrix[n])
378 ofs << i << '\t' << j << 333 {
379 << 334 ofs << G4endl;
380 // Total dose << 335 ofs << i << '\t' << j << '\t' << k << '\t';
381 //ofs << std::setw(wid << 336 // Total dose
382 ofs << (matrix[n]/mass << 337 ofs << std::setw(width) << matrix[n]/massOfVoxel/doseUnit;
383 << 338 if (secondary)
384 if (secondary) << 339 {
385 { << 340 for (size_t l=0; l < ionStore.size(); l++)
386 for (size_t l=0; l << 341 {
387 { << 342 // Fill ASCII file rows
388 // Fill ASCII << 343 ofs << std::setw(width) << ionStore[l].dose[n]/massOfVoxel/doseUnit <<
389 //ofs << std:: << 344 std::setw(width) << ionStore[l].fluence[n];
390 // std:: << 345 }
391 << 346 }
392 ofs << '\t' < << 347 }
393 '\t' << i << 348 }
394 } << 349 ofs.close();
395 } << 350 }
396 } << 351 }
397 } << 352 /////////////////////////////////////////////////////////////////////////////
398 ofs.close(); << 353
399 } << 354 #ifdef G4ANALYSIS_USE_ROOT
400 } << 355 void HadrontherapyMatrix::StoreDoseFluenceRoot()
401 ////////////////////////////////////////////// <<
402 void HadrontherapyMatrix::Fill(G4int i, G4int <<
403 G4double energy <<
404 { 356 {
405 if (matrix) << 357 HadrontherapyAnalysisManager* analysis = HadrontherapyAnalysisManager::GetInstance();
406 matrix[Index(i,j,k)] += energyDeposit; << 358 if (analysis -> IsTheTFile())
>> 359 {
>> 360 for(G4int i = 0; i < numberOfVoxelAlongX; i++)
>> 361 for(G4int j = 0; j < numberOfVoxelAlongY; j++)
>> 362 for(G4int k = 0; k < numberOfVoxelAlongZ; k++)
>> 363 {
>> 364 G4int n = Index(i, j, k);
>> 365 for (size_t l=0; l < ionStore.size(); l++)
>> 366
>> 367 {
>> 368 // Do the same work for .root file: fill dose/fluence ntuple
>> 369 analysis -> FillVoxelFragmentTuple( i, j, k,
>> 370 ionStore[l].A,
>> 371 ionStore[l].Z,
>> 372 ionStore[l].dose[n]/massOfVoxel/doseUnit,
>> 373 ionStore[l].fluence[n] );
407 374
408 // Store the energy deposit in the matrix <<
409 // to the phantom voxel <<
410 } <<
411 375
>> 376 }
>> 377 }
>> 378 }
>> 379 }
>> 380 #endif
412 381
>> 382 void HadrontherapyMatrix::Fill(G4int i, G4int j, G4int k,
>> 383 G4double energyDeposit)
>> 384 {
>> 385 if (matrix)
>> 386 matrix[Index(i,j,k)] += energyDeposit;
>> 387
>> 388 // Store the energy deposit in the matrix element corresponding
>> 389 // to the phantom voxel
>> 390 }
>> 391 void HadrontherapyMatrix::TotalEnergyDeposit()
>> 392 {
>> 393 // Convert energy deposited to dose.
>> 394 // Store the information of the matrix in a ntuple and in
>> 395 // a 1D Histogram
>> 396 #ifdef G4ANALYSIS_USE_ROOT
>> 397 HadrontherapyAnalysisManager* analysis = HadrontherapyAnalysisManager::GetInstance();
>> 398 #endif
413 399
>> 400 if (matrix)
>> 401 {
>> 402 for(G4int i = 0; i < numberOfVoxelAlongX; i++)
>> 403 for(G4int j = 0; j < numberOfVoxelAlongY; j++)
>> 404 for(G4int k = 0; k < numberOfVoxelAlongZ; k++)
>> 405 {
>> 406 #ifdef G4ANALYSIS_USE_ROOT
>> 407 G4int n = Index(i,j,k);
>> 408 if (analysis -> IsTheTFile() )
>> 409 {
>> 410 analysis -> FillEnergyDeposit(i, j, k, matrix[n]/massOfVoxel/doseUnit);
>> 411 analysis -> BraggPeak(i, matrix[n]/massOfVoxel/doseUnit);
>> 412 }
>> 413 #endif
>> 414 }
>> 415 }
>> 416 }
414 417
415 418