Warning, file /geant4/examples/advanced/hadrontherapy/src/HadrontherapyMatrix.cc was not indexed
or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029 #include <fstream>
0030 #include <iostream>
0031 #include <sstream>
0032 #include <iomanip>
0033
0034 #include "HadrontherapyMatrix.hh"
0035 #include "HadrontherapyPrimaryGeneratorAction.hh"
0036 #include "globals.hh"
0037 #include "G4SystemOfUnits.hh"
0038 #include "G4RunManager.hh"
0039 #include "G4ParticleGun.hh"
0040 #include "HadrontherapySteppingAction.hh"
0041 #include "HadrontherapyAnalysisFileMessenger.hh"
0042 #include "G4SystemOfUnits.hh"
0043 #include <time.h>
0044
0045 HadrontherapyAnalysis* HadrontherapyAnalysis::instance = 0;
0046
0047
0048 HadrontherapyAnalysis::HadrontherapyAnalysis()
0049 {
0050 fMess = new HadrontherapyAnalysisFileMessenger(this);
0051 }
0052
0053
0054 HadrontherapyAnalysis::~HadrontherapyAnalysis()
0055 {
0056 delete fMess;
0057 }
0058
0059
0060 HadrontherapyAnalysis* HadrontherapyAnalysis::GetInstance(){
0061
0062 if (instance == 0) instance = new HadrontherapyAnalysis;
0063 return instance;
0064 }
0065
0066 HadrontherapyMatrix* HadrontherapyMatrix::instance = NULL;
0067 G4bool HadrontherapyMatrix::secondary = false;
0068
0069
0070
0071 HadrontherapyMatrix* HadrontherapyMatrix::GetInstance()
0072 {
0073 return instance;
0074 }
0075
0076
0077
0078
0079 HadrontherapyMatrix* HadrontherapyMatrix::GetInstance(G4int voxelX, G4int voxelY, G4int voxelZ, G4double mass)
0080 {
0081 if (instance) delete instance;
0082 instance = new HadrontherapyMatrix(voxelX, voxelY, voxelZ, mass);
0083 instance -> Initialize();
0084 return instance;
0085 }
0086
0087
0088 HadrontherapyMatrix::HadrontherapyMatrix(G4int voxelX, G4int voxelY, G4int voxelZ, G4double mass):
0089 stdFile("Dose.out"),
0090 doseUnit(gray)
0091 {
0092
0093
0094
0095 numberOfVoxelAlongX = voxelX;
0096 numberOfVoxelAlongY = voxelY;
0097 numberOfVoxelAlongZ = voxelZ;
0098 massOfVoxel = mass;
0099
0100
0101
0102 matrix = new G4double[numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ];
0103 if (matrix)
0104 {
0105 G4cout << "HadrontherapyMatrix: Memory space to store physical dose into " <<
0106 numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ <<
0107 " voxels has been allocated " << G4endl;
0108 }
0109
0110 else G4Exception("HadrontherapyMatrix::HadrontherapyMatrix()", "Hadrontherapy0005", FatalException, "Can't allocate memory to store physical dose!");
0111
0112
0113
0114
0115
0116
0117 hitTrack = new G4int[numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ];
0118 ClearHitTrack();
0119 }
0120
0121
0122 HadrontherapyMatrix::~HadrontherapyMatrix()
0123 {
0124 delete[] matrix;
0125 delete[] hitTrack;
0126 Clear();
0127 }
0128
0129
0130 void HadrontherapyMatrix::Clear()
0131 {
0132 for (size_t i=0; i<ionStore.size(); i++)
0133 {
0134 delete[] ionStore[i].dose;
0135 delete[] ionStore[i].fluence;
0136 }
0137 ionStore.clear();
0138 }
0139
0140
0141
0142
0143 void HadrontherapyMatrix::Initialize()
0144 {
0145
0146 Clear();
0147
0148 for(int i=0;i<numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ;i++)
0149 {
0150 matrix[i] = 0;
0151 }
0152 }
0153
0154
0155
0156
0157
0158 void HadrontherapyMatrix::PrintNuclides()
0159 {
0160 for (size_t i=0; i<ionStore.size(); i++)
0161 {
0162 G4cout << ionStore[i].name << G4endl;
0163 }
0164 }
0165
0166
0167
0168
0169 void HadrontherapyMatrix::ClearHitTrack()
0170 {
0171 for(G4int i=0; i<numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ; i++) hitTrack[i] = 0;
0172 }
0173
0174
0175 G4int* HadrontherapyMatrix::GetHitTrack(G4int i, G4int j, G4int k)
0176 {
0177 return &(hitTrack[Index(i,j,k)]);
0178 }
0179
0180
0181
0182
0183
0184
0185
0186 G4bool HadrontherapyMatrix::Fill(G4int trackID,
0187 G4ParticleDefinition* particleDef,
0188 G4int i, G4int j, G4int k,
0189 G4double energyDeposit,
0190 G4bool fluence)
0191 {
0192
0193 if ( (energyDeposit <=0. && !fluence) || !secondary) return false;
0194
0195
0196 G4int PDGencoding = particleDef -> GetPDGEncoding();
0197 PDGencoding -= PDGencoding%10;
0198
0199
0200 for (size_t l=0; l < ionStore.size(); l++)
0201 {
0202 if (ionStore[l].PDGencoding == PDGencoding )
0203 {
0204
0205 if ( (trackID ==1 && ionStore[l].isPrimary) || (trackID !=1 && !ionStore[l].isPrimary))
0206 {
0207 if (energyDeposit > 0.)
0208
0209 ionStore[l].dose[Index(i, j, k)] += energyDeposit;
0210
0211
0212
0213 if (fluence) ionStore[l].fluence[Index(i, j, k)]++;
0214 return true;
0215 }
0216 }
0217 }
0218 G4int Z = particleDef-> GetAtomicNumber();
0219 G4int A = particleDef-> GetAtomicMass();
0220 G4String fullName = particleDef -> GetParticleName();
0221 G4String name = fullName.substr (0, fullName.find("[") );
0222
0223
0224 ion newIon =
0225 {
0226 (trackID == 1) ? true:false,
0227 PDGencoding,
0228 name,
0229 name.length(),
0230 Z,
0231 A,
0232 new G4double[numberOfVoxelAlongX * numberOfVoxelAlongY * numberOfVoxelAlongZ],
0233 new unsigned int[numberOfVoxelAlongX * numberOfVoxelAlongY * numberOfVoxelAlongZ]
0234 };
0235
0236
0237
0238 if (newIon.dose && newIon.fluence)
0239 {
0240 for(G4int q=0; q<numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ; q++)
0241 {
0242 newIon.dose[q] = 0.;
0243 newIon.fluence[q] = 0;
0244 }
0245
0246 if (energyDeposit > 0.) newIon.dose[Index(i, j, k)] += energyDeposit;
0247 if (fluence) newIon.fluence[Index(i, j, k)]++;
0248
0249 ionStore.push_back(newIon);
0250 return true;
0251 }
0252
0253 else
0254
0255 {
0256 return false;
0257 }
0258 }
0259
0260
0261
0262
0263
0264
0265
0266
0267 void HadrontherapyMatrix::StoreMatrix(G4String file, void* data, size_t psize)
0268 {
0269 if (data)
0270 {
0271 ofs.open(file, std::ios::out);
0272 if (ofs.is_open())
0273 {
0274 for(G4int i = 0; i < numberOfVoxelAlongX; i++)
0275 for(G4int j = 0; j < numberOfVoxelAlongY; j++)
0276 for(G4int k = 0; k < numberOfVoxelAlongZ; k++)
0277 {
0278 G4int n = Index(i, j, k);
0279
0280 if (psize == sizeof(unsigned int))
0281 {
0282 unsigned int* pdata = (unsigned int*)data;
0283
0284 if (pdata[n])
0285
0286 ofs << i << '\t' << j << '\t' << k << '\t' << pdata[n] << G4endl;
0287 }
0288
0289 else if (psize == sizeof(G4double))
0290
0291 {
0292 G4double* pdata = (G4double*)data;
0293 if (pdata[n]) ofs << i << '\t' << j << '\t' << k << '\t' << pdata[n] << G4endl;
0294 }
0295 }
0296 ofs.close();
0297 }
0298 }
0299 }
0300
0301
0302
0303 void HadrontherapyMatrix::StoreFluenceData()
0304 {
0305 for (size_t i=0; i < ionStore.size(); i++){
0306 StoreMatrix(ionStore[i].name + "_Fluence.out", ionStore[i].fluence, sizeof(unsigned int));
0307 }
0308 }
0309
0310
0311
0312 void HadrontherapyMatrix::StoreDoseData()
0313 {
0314
0315 for (size_t i=0; i < ionStore.size(); i++){
0316 StoreMatrix(ionStore[i].name + "_Dose.out", ionStore[i].dose, sizeof(G4double));
0317 }
0318 }
0319
0320
0321
0322
0323
0324
0325 void HadrontherapyMatrix::StoreDoseFluenceAscii(G4String file)
0326 {
0327 #define width 15L
0328 filename = (file=="") ? stdFile:file;
0329
0330
0331
0332 std::sort(ionStore.begin(), ionStore.end());
0333 G4cout << "Dose is being written to " << filename << G4endl;
0334 ofs.open(filename, std::ios::out);
0335
0336 if (ofs.is_open())
0337 {
0338
0339
0340 ofs << "i" << '\t' << "j" << '\t' << "k";
0341 G4cout << "i" << '\t' << "j" << '\t' << "k";
0342
0343
0344 ofs <<'\t' <<"Dose(Gy)";
0345
0346 G4cout << '\t' << "Dose(Gy)";
0347
0348 G4String fluence = "_f";
0349 if (secondary)
0350 {
0351 for (size_t l=0; l < ionStore.size(); l++)
0352 {
0353 G4String a = (ionStore[l].isPrimary) ? "_1":"";
0354
0355
0356
0357
0358 ofs << '\t' << ionStore[l].name + a <<
0359 '\t' << ionStore[l].name + a + fluence;
0360
0361 G4cout << '\t' << ionStore[l].name + a <<
0362 '\t' << ionStore[l].name + a + fluence;
0363
0364 }
0365
0366 }
0367
0368
0369 for(G4int i = 0; i < numberOfVoxelAlongX; i++)
0370 for(G4int j = 0; j < numberOfVoxelAlongY; j++)
0371 for(G4int k = 0; k < numberOfVoxelAlongZ; k++)
0372 {
0373 G4int n = Index(i, j, k);
0374
0375 if (matrix[n])
0376 {
0377 ofs << G4endl;
0378 ofs << i << '\t' << j << '\t' << k << '\t';
0379
0380
0381
0382 ofs << (matrix[n]/massOfVoxel)/doseUnit;
0383
0384 if (secondary)
0385 {
0386 for (size_t l=0; l < ionStore.size(); l++)
0387 {
0388
0389
0390
0391
0392 ofs << '\t' << ionStore[l].dose[n]/massOfVoxel/doseUnit <<
0393 '\t' << ionStore[l].fluence[n];
0394 }
0395 }
0396 }
0397 }
0398 ofs.close();
0399 }
0400 }
0401
0402 void HadrontherapyMatrix::Fill(G4int i, G4int j, G4int k,
0403 G4double energyDeposit)
0404 {
0405 if (matrix)
0406 matrix[Index(i,j,k)] += energyDeposit;
0407
0408
0409
0410 }
0411
0412
0413
0414