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