File indexing completed on 2025-02-26 09:19:38
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