- Timestamp:
- Jun 14, 2010, 3:54:58 PM (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/examples/advanced/hadrontherapy/src/HadrontherapyMatrix.cc
r1230 r1313 1 //2 // ********************************************************************3 // * License and Disclaimer *4 // * *5 // * The Geant4 software is copyright of the Copyright Holders of *6 // * the Geant4 Collaboration. It is provided under the terms and *7 // * conditions of the Geant4 Software License, included in the file*8 // * LICENSE and available at http://cern.ch/geant4/license . These *9 // * include a list of copyright holders. *10 // * *11 // * Neither the authors of this software system, nor their employing *12 // * institutes,nor the agencies providing financial support for this *13 // * work make any representation or warranty, express or implied, *14 // * regarding this software system or assume any liability for its *15 // * use. Please see the license in the file LICENSE and URL above *16 // * for the full disclaimer and the limitation of liability. *17 // * *18 // * This code implementation is the result of the scientific and *19 // * technical work of the GEANT4 collaboration. *20 // * By using, copying, modifying or distributing the software (or *21 // * any work based on the software) you agree to acknowledge its *22 // * use in resulting scientific publications, and indicate your *23 // * acceptance of all terms of the Geant4 Software license. *24 // ********************************************************************25 //26 // $Id: HadrontherapyMatrix.cc;27 // See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy//1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file* 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 // $Id: HadrontherapyMatrix.cc; 27 // See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy// 28 28 29 29 #include "HadrontherapyMatrix.hh" 30 30 #include "HadrontherapyAnalysisManager.hh" 31 #include "G4RunManager.hh" 32 #include "HadrontherapyPrimaryGeneratorAction.hh" 33 #include "G4ParticleGun.hh" 34 35 #include <fstream> 36 #include <unistd.h> 37 #include <iostream> 38 #include <sstream> 39 #include <iomanip> 31 40 #include "globals.hh" 32 #include <fstream> 33 41 42 // Units definition: CLHEP/Units/SystemOfUnits.h 43 // 34 44 HadrontherapyMatrix* HadrontherapyMatrix::instance = NULL; 35 // Static method that only return a pointer to the matrix object 36 HadrontherapyMatrix* HadrontherapyMatrix::getInstance() 37 { 38 return instance; 39 } 40 // This STATIC method delete (!) the old matrix and rewrite a new object returning a pointer to it 41 HadrontherapyMatrix* HadrontherapyMatrix::getInstance(G4int voxelX, G4int voxelY, G4int voxelZ) 42 { 43 if (instance) delete instance; 44 instance = new HadrontherapyMatrix(voxelX, voxelY, voxelZ); 45 instance -> Initialize(); 46 return instance; 47 } 48 49 HadrontherapyMatrix::HadrontherapyMatrix(G4int voxelX, G4int voxelY, G4int voxelZ): 50 matrix(0) 45 G4bool HadrontherapyMatrix::secondaries = false; 46 // Only return a pointer to matrix 47 HadrontherapyMatrix* HadrontherapyMatrix::GetInstance() 48 { 49 return instance; 50 } 51 // This STATIC method delete (!) the old matrix and rewrite a new object returning a pointer to it 52 // TODO A check on the parameters is required! 53 HadrontherapyMatrix* HadrontherapyMatrix::GetInstance(G4int voxelX, G4int voxelY, G4int voxelZ, G4double mass) 54 { 55 if (instance) delete instance; 56 instance = new HadrontherapyMatrix(voxelX, voxelY, voxelZ, mass); 57 instance -> Initialize();// XXX 58 return instance; 59 } 60 HadrontherapyMatrix::HadrontherapyMatrix(G4int voxelX, G4int voxelY, G4int voxelZ, G4double mass): 61 filename("Dose.out"), 62 doseUnit(MeV/g) 51 63 { 52 // Number of the voxels of the phantom 53 // For Y = Z = 1 the phantom is divided in slices (and not in voxels) 54 // orthogonal to the beam axis 55 numberVoxelX = voxelX; 56 numberVoxelY = voxelY; 57 numberVoxelZ = voxelZ; 58 // Create the dose matrix 59 matrix = new G4double[numberVoxelX*numberVoxelY*numberVoxelZ]; 60 if (matrix) G4cout << "Matrix: Memory space to store physical dose into " << 61 numberVoxelX*numberVoxelY*numberVoxelZ << 62 " voxels has been allocated " << G4endl; 63 else G4Exception("Can't allocate memory to store physical dose!"); 64 } 65 64 // Number of the voxels of the phantom 65 // For Y = Z = 1 the phantom is divided in slices (and not in voxels) 66 // orthogonal to the beam axis 67 numberOfVoxelAlongX = voxelX; 68 numberOfVoxelAlongY = voxelY; 69 numberOfVoxelAlongZ = voxelZ; 70 massOfVoxel = mass; 71 // Create the dose matrix 72 matrix = new G4double[numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ]; 73 if (matrix) 74 { 75 G4cout << "HadrontherapyMatrix: Memory space to store physical dose into " << 76 numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ << 77 " voxels has been allocated " << G4endl; 78 } 79 else G4Exception(" HadrontherapyMatrix::HadrontherapyMatrix. Can't allocate memory to store physical dose!"); 80 // Hit voxel (TrackID) marker 81 // This array mark the status of voxel, if a hit occur, with the trackID of the particle 82 // Must be initialized 83 hitTrack = new G4int[numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ]; 84 ClearHitTrack(); 85 } 86 87 ///////////////////////////////////////////////////////////////////////////// 66 88 HadrontherapyMatrix::~HadrontherapyMatrix() 67 89 { 68 90 delete[] matrix; 69 } 70 71 void HadrontherapyMatrix::flush() 72 { 73 if(matrix) 74 for(int i=0; i<numberVoxelX*numberVoxelY*numberVoxelZ; i++) 75 { 76 matrix[i] = 0; 77 } 78 } 91 delete[] hitTrack; 92 // clear fluences/dose data 93 Clear(); 94 } 95 96 ///////////////////////////////////////////////////////////////////////////// 97 void HadrontherapyMatrix::Clear() 98 { 99 for (size_t i=0; i<ionStore.size(); i++) 100 { 101 delete[] ionStore[i].dose; 102 delete[] ionStore[i].fluence; 103 } 104 ionStore.clear(); 105 } 106 107 ///////////////////////////////////////////////////////////////////////////// 108 // Initialise the elements of the matrix to zero 79 109 void HadrontherapyMatrix::Initialize() 80 110 { 81 // Initialise the elements of the matrix to zero 82 for(G4int i = 0; i < numberVoxelX; i++) 83 { 84 for(G4int j = 0; j < numberVoxelY; j++) 85 { 86 for(G4int k = 0; k < numberVoxelZ; k++) 87 88 matrix[Index(i,j,k)] = 0.; 89 } 90 } 91 } 111 Clear(); 112 for(int i=0;i<numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ;i++) 113 { 114 matrix[i] = 0; 115 } 116 } 117 ///////////////////////////////////////////////////////////////////////////// 118 ///////////////////////////////////////////////////////////////////////////// 119 // Print generated nuclides list 120 void HadrontherapyMatrix::PrintNuclides() 121 { 122 for (size_t i=0; i<ionStore.size(); i++) 123 { 124 G4cout << ionStore[i].name << G4endl; 125 } 126 } 127 ///////////////////////////////////////////////////////////////////////////// 128 // Clear Hit voxel (TrackID) markers 129 void HadrontherapyMatrix::ClearHitTrack() 130 { 131 for(G4int i=0; i<numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ; i++) hitTrack[i] = 0; 132 } 133 // Return Hit status 134 G4int* HadrontherapyMatrix::GetHitTrack(G4int i, G4int j, G4int k) 135 { 136 return &(hitTrack[Index(i,j,k)]); 137 } 138 ///////////////////////////////////////////////////////////////////////////// 139 // Dose methods... 140 // Fill DOSE/fluence matrix for particle/ion: 141 // If fluence parameter is true (default value is FALSE) then fluence at voxel (i, j, k) is increased. 142 // The energyDeposit parameter fill the dose matrix for voxel (i,j,k) 143 ///////////////////////////////////////////////////////////////////////////// 144 145 G4bool HadrontherapyMatrix::Fill(G4int trackID, 146 G4ParticleDefinition* particleDef, 147 G4int i, G4int j, G4int k, 148 G4double energyDeposit, 149 G4bool fluence) 150 { 151 if (energyDeposit <=0. && !fluence || !secondaries) return false; 152 // Get Particle Data Group particle ID 153 G4int PDGencoding = particleDef -> GetPDGEncoding(); 154 PDGencoding -= PDGencoding%10; 155 156 // Search for already allocated data... 157 for (size_t l=0; l < ionStore.size(); l++) 158 { 159 if (ionStore[l].PDGencoding == PDGencoding ) 160 { // Is it a primary or a secondary particle? 161 if ( trackID ==1 && ionStore[l].isPrimary || trackID !=1 && !ionStore[l].isPrimary) 162 { 163 if (energyDeposit > 0.) ionStore[l].dose[Index(i, j, k)] += energyDeposit/massOfVoxel; 164 165 // Fill a matrix per each ion with the fluence 166 if (fluence) ionStore[l].fluence[Index(i, j, k)]++; 167 return true; 168 } 169 } 170 } 171 172 G4int Z = particleDef-> GetAtomicNumber(); 173 G4int A = particleDef-> GetAtomicMass(); 174 175 G4String fullName = particleDef -> GetParticleName(); 176 G4String name = fullName.substr (0, fullName.find("[") ); // cut excitation energy 177 // Let's put a new particle in our store... 178 ion newIon = 179 { 180 (trackID == 1) ? true:false, 181 PDGencoding, 182 name, 183 name.length(), 184 Z, 185 A, 186 new G4double[numberOfVoxelAlongX * numberOfVoxelAlongY * numberOfVoxelAlongZ], 187 new unsigned int[numberOfVoxelAlongX * numberOfVoxelAlongY * numberOfVoxelAlongZ] 188 }; 189 // Initialize data 190 if (newIon.dose && newIon.fluence) 191 { 192 for(G4int m=0; m<numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ; m++) 193 { 194 newIon.dose[m] = 0.; 195 newIon.fluence[m] = 0; 196 } 197 if (energyDeposit > 0.) newIon.dose[Index(i, j, k)] += energyDeposit/massOfVoxel; 198 if (fluence) newIon.fluence[Index(i, j, k)]++; 199 200 ionStore.push_back(newIon); 201 202 // TODO Put some verbosity check 203 /* 204 G4cout << "Memory space to store the DOSE/FLUENCE into " << 205 numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ << 206 " voxels has been allocated for the nuclide " << newIon.name << 207 " (Z = " << Z << ", A = " << A << ")" << G4endl ; 208 */ 209 return true; 210 } 211 else // XXX Out of memory! XXX 212 { 213 return false; 214 } 215 216 } 217 218 ///////////////////////////////////////////////////////////////////////////// 219 ///////////////////////////////////////////////////////////////////////////// 220 // Methods to store data to filenames... 221 //////////////////////////////////////////////////////////////////////////// 222 //////////////////////////////////////////////////////////////////////////// 223 // 224 // General method to store matrix data to filename 225 void HadrontherapyMatrix::StoreMatrix(G4String file, void* data, size_t psize) 226 { 227 if (data) 228 { 229 ofs.open(file, std::ios::out); 230 if (ofs.is_open()) 231 { 232 for(G4int i = 0; i < numberOfVoxelAlongX; i++) 233 for(G4int j = 0; j < numberOfVoxelAlongY; j++) 234 for(G4int k = 0; k < numberOfVoxelAlongZ; k++) 235 { 236 G4int n = Index(i, j, k); 237 // Check for data type: u_int, G4double, XXX 238 if (psize == sizeof(unsigned int)) 239 { 240 unsigned int* pdata = (unsigned int*)data; 241 if (pdata[n]) ofs << i << '\t' << j << '\t' << 242 k << '\t' << pdata[n] << G4endl; 243 } 244 else if (psize == sizeof(G4double)) 245 { 246 G4double* pdata = (G4double*)data; 247 if (pdata[n]) ofs << i << '\t' << j << '\t' << 248 k << '\t' << pdata[n] << G4endl; 249 } 250 } 251 ofs.close(); 252 } 253 } 254 } 255 256 // Store fluence per single ion in multiple files 257 void HadrontherapyMatrix::StoreFluenceData() 258 { 259 for (size_t i=0; i < ionStore.size(); i++){ 260 StoreMatrix(ionStore[i].name + "_Fluence.out", ionStore[i].fluence, sizeof(unsigned int)); 261 } 262 } 263 // Store dose per single ion in multiple files 264 void HadrontherapyMatrix::StoreDoseData() 265 { 266 267 for (size_t i=0; i < ionStore.size(); i++){ 268 StoreMatrix(ionStore[i].name + "_Dose.out", ionStore[i].dose, sizeof(G4double)); 269 } 270 } 271 ///////////////////////////////////////////////////////////////////////// 272 // Store dose for all ions into a single file and into ntuples. 273 // Please note that this function is called via messenger commands 274 // defined in the HadrontherapyAnalysisFileMessenger.cc class file 275 void HadrontherapyMatrix::StoreDoseFluenceAscii() 276 { 277 // Sort like periodic table 278 std::sort(ionStore.begin(), ionStore.end()); 279 #define width 15L 280 ofs.open(filename, std::ios::out); 281 if (ofs.is_open()) 282 { 283 284 // 285 // Write the voxels index and the list of particles/ions 286 ofs << std::setprecision(6) << std::left << 287 "i\tj\tk\t"; 288 /* 289 G4RunManager *runManager = G4RunManager::GetRunManager(); 290 HadrontherapyPrimaryGeneratorAction *pPGA = (HadrontherapyPrimaryGeneratorAction*)runManager -> GetUserPrimaryGeneratorAction(); 291 G4String name = pPGA -> GetParticleGun() -> GetParticleDefinition() -> GetParticleName(); 292 */ 293 // Total dose 294 ofs << std::setw(width) << "Dose"; 295 296 for (size_t l=0; l < ionStore.size(); l++) 297 { 298 G4String a = (ionStore[l].isPrimary) ? "_1":""; // is it a primary? 299 ofs << std::setw(width) << ionStore[l].name + a << 300 std::setw(width) << ionStore[l].name + a; 301 } 302 ofs << G4endl; 303 304 /* 305 // PDGencoding 306 ofs << std::setprecision(6) << std::left << 307 "0\t0\t0\t"; 308 309 // Total dose 310 ofs << std::setw(width) << '0'; 311 for (size_t l=0; l < ionStore.size(); l++) 312 { 313 ofs << std::setw(width) << ionStore[l].PDGencoding << 314 std::setw(width) << ionStore[l].PDGencoding; 315 } 316 ofs << G4endl; 317 */ 318 // Write data 319 for(G4int i = 0; i < numberOfVoxelAlongX; i++) 320 for(G4int j = 0; j < numberOfVoxelAlongY; j++) 321 for(G4int k = 0; k < numberOfVoxelAlongZ; k++) 322 { 323 G4int n = Index(i, j, k); 324 // Write only not identically null data lines 325 if (matrix[n]) 326 { 327 ofs << G4endl; 328 ofs << i << '\t' << j << '\t' << k << '\t'; 329 // Total dose 330 ofs << std::setw(width) << matrix[n]/(doseUnit); 331 { 332 for (size_t l=0; l < ionStore.size(); l++) 333 { 334 // Fill ASCII file rows 335 ofs << std::setw(width) << ionStore[l].dose[n]/(doseUnit) << 336 std::setw(width) << ionStore[l].fluence[n]; 337 } 338 } 339 } 340 } 341 ofs.close(); 342 } 343 } 344 ///////////////////////////////////////////////////////////////////////////// 345 346 #ifdef G4ANALYSIS_USE_ROOT 347 void HadrontherapyMatrix::StoreDoseFluenceRoot() 348 { 349 HadrontherapyAnalysisManager* analysis = HadrontherapyAnalysisManager::GetInstance(); 350 for(G4int i = 0; i < numberOfVoxelAlongX; i++) 351 for(G4int j = 0; j < numberOfVoxelAlongY; j++) 352 for(G4int k = 0; k < numberOfVoxelAlongZ; k++) 353 { 354 G4int n = Index(i, j, k); 355 for (size_t l=0; l < ionStore.size(); l++) 356 357 { 358 // Do the same work for .root file: fill dose/fluence ntuple 359 analysis -> FillVoxelFragmentTuple( i, j, k, 360 ionStore[l].A, 361 ionStore[l].Z, 362 ionStore[l].dose[n]/(doseUnit), 363 ionStore[l].fluence[n] ); 364 365 366 } 367 } 368 } 369 #endif 92 370 93 371 void HadrontherapyMatrix::Fill(G4int i, G4int j, G4int k, 94 372 G4double energyDeposit) 95 373 { 96 if (matrix) 97 matrix[Index(i,j,k)] += energyDeposit; 98 99 // Store the energy deposit in the matrix element corresponding 100 // to the phantom voxel i, j, k 101 } 102 374 if (matrix) 375 matrix[Index(i,j,k)] += energyDeposit; 376 377 // Store the energy deposit in the matrix element corresponding 378 // to the phantom voxel 379 } 103 380 void HadrontherapyMatrix::TotalEnergyDeposit() 104 381 { 382 // Convert energy deposited to dose. 105 383 // Store the information of the matrix in a ntuple and in 106 384 // a 1D Histogram 107 385 108 G4int k; 109 G4int j; 110 G4int i; 111 112 if (matrix) 386 #ifdef G4ANALYSIS_USE_ROOT 387 HadrontherapyAnalysisManager* analysis = HadrontherapyAnalysisManager::GetInstance(); 388 #endif 389 if (matrix) 113 390 { 114 std::ofstream ofs; 115 ofs.open("DoseDistribution.out"); 116 117 for(G4int l = 0; l < numberVoxelZ; l++) 118 { 119 k = l; 120 121 for(G4int m = 0; m < numberVoxelY; m++) 122 { 123 j = m * numberVoxelZ + k; 124 125 for(G4int n = 0; n < numberVoxelX; n++) 126 { 127 i = n* numberVoxelZ * numberVoxelY + j; 128 if(matrix[i] != 0) 129 { 130 ofs << n << '\t' << m << '\t' << 131 k << '\t' << matrix[i] << G4endl; 132 133 #ifdef ANALYSIS_USE 134 HadrontherapyAnalysisManager* analysis = 135 HadrontherapyAnalysisManager::getInstance(); 136 analysis -> FillEnergyDeposit(n, m, k, matrix[i]); 137 analysis -> BraggPeak(n, matrix[i]); 391 for(G4int i = 0; i < numberOfVoxelAlongX; i++) 392 for(G4int j = 0; j < numberOfVoxelAlongY; j++) 393 for(G4int k = 0; k < numberOfVoxelAlongZ; k++) 394 { 395 G4int n = Index(i,j,k); 396 matrix[n] = matrix[n] / massOfVoxel; 397 #ifdef G4ANALYSIS_USE_ROOT 398 analysis -> FillEnergyDeposit(i, j, k, matrix[n]/doseUnit); 399 analysis -> BraggPeak(i, matrix[n]/doseUnit); 138 400 #endif 139 } 140 } 141 } 142 } 143 } 144 } 401 } 402 } 403 } 404
Note: See TracChangeset
for help on using the changeset viewer.