Geant4 Cross Reference |
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 // Hadrontherapy advanced example for Geant4 27 // See more at: https://twiki.cern.ch/twiki/bin/view/Geant4/AdvancedExamplesHadrontherapy 28 29 #include "HadrontherapyDetectorConstruction.hh" 30 #include "HadrontherapyLet.hh" 31 32 #include "HadrontherapyMatrix.hh" 33 #include "HadrontherapyInteractionParameters.hh" 34 #include "HadrontherapyPrimaryGeneratorAction.hh" 35 #include "HadrontherapyMatrix.hh" 36 #include "G4AnalysisManager.hh" 37 #include "G4RunManager.hh" 38 #include "G4SystemOfUnits.hh" 39 #include <cmath> 40 41 HadrontherapyLet* HadrontherapyLet::instance = NULL; 42 G4bool HadrontherapyLet::doCalculation = false; 43 44 HadrontherapyLet* HadrontherapyLet::GetInstance(HadrontherapyDetectorConstruction *pDet) 45 { 46 if (instance) delete instance; 47 instance = new HadrontherapyLet(pDet); 48 return instance; 49 } 50 51 HadrontherapyLet* HadrontherapyLet::GetInstance() 52 { 53 return instance; 54 } 55 56 HadrontherapyLet::HadrontherapyLet(HadrontherapyDetectorConstruction* pDet) 57 :filename("Let.out"),matrix(0) // Default output filename 58 { 59 60 matrix = HadrontherapyMatrix::GetInstance(); 61 62 if (!matrix) 63 G4Exception("HadrontherapyLet::HadrontherapyLet", 64 "Hadrontherapy0005", FatalException, 65 "HadrontherapyMatrix not found. Firstly create an instance of it."); 66 67 nVoxels = matrix -> GetNvoxel(); 68 69 numberOfVoxelAlongX = matrix -> GetNumberOfVoxelAlongX(); 70 numberOfVoxelAlongY = matrix -> GetNumberOfVoxelAlongY(); 71 numberOfVoxelAlongZ = matrix -> GetNumberOfVoxelAlongZ(); 72 73 G4RunManager *runManager = G4RunManager::GetRunManager(); 74 pPGA = (HadrontherapyPrimaryGeneratorAction*)runManager -> GetUserPrimaryGeneratorAction(); 75 // Pointer to the detector material 76 detectorMat = pDet -> GetDetectorLogicalVolume() -> GetMaterial(); 77 density = detectorMat -> GetDensity(); 78 // Instances for Total LET 79 totalLetD = new G4double[nVoxels]; 80 DtotalLetD = new G4double[nVoxels]; 81 totalLetT = new G4double[nVoxels]; 82 DtotalLetT = new G4double[nVoxels]; 83 84 } 85 86 HadrontherapyLet::~HadrontherapyLet() 87 { 88 Clear(); 89 delete [] totalLetD; 90 delete [] DtotalLetD; 91 delete [] totalLetT; 92 delete [] DtotalLetT; 93 } 94 95 // Fill energy spectrum for every voxel (local energy spectrum) 96 void HadrontherapyLet::Initialize() 97 { 98 for(G4int v=0; v < nVoxels; v++) totalLetD[v] = DtotalLetD[v] = totalLetT[v] = DtotalLetT[v] = 0.; 99 Clear(); 100 } 101 /** 102 * Clear all stored data 103 */ 104 void HadrontherapyLet::Clear() 105 { 106 for (size_t i=0; i < ionLetStore.size(); i++) 107 { 108 delete [] ionLetStore[i].letDN; // Let Dose Numerator 109 delete [] ionLetStore[i].letDD; // Let Dose Denominator 110 delete [] ionLetStore[i].letTN; // Let Track Numerator 111 delete [] ionLetStore[i].letTD; // Let Track Denominator 112 } 113 ionLetStore.clear(); 114 } 115 void HadrontherapyLet::FillEnergySpectrum(G4int trackID, 116 G4ParticleDefinition* particleDef, 117 G4double ekinMean, 118 const G4Material* mat, 119 G4double DE, 120 G4double DEEletrons, 121 G4double DX, 122 G4int i, G4int j, G4int k) 123 { 124 if (DE <= 0. || DX <=0. ) return; // calculate only energy deposit 125 if (!doCalculation) return; 126 127 // atomic number 128 G4int Z = particleDef -> GetAtomicNumber(); 129 if (Z<1) return; // calculate only protons and ions 130 131 G4int PDGencoding = particleDef -> GetPDGEncoding(); 132 PDGencoding -= PDGencoding%10; // simple Particle data group id without excitation level 133 134 G4int voxel = matrix -> Index(i,j,k); 135 136 // ICRU stopping power calculation 137 G4EmCalculator emCal; 138 // use the mean kinetic energy of ions in a step to calculate ICRU stopping power 139 G4double Lsn = emCal.ComputeElectronicDEDX(ekinMean, particleDef, mat); 140 141 142 // Total LET calculation... 143 totalLetD[voxel] += (DE + DEEletrons) * Lsn; // total dose Let Numerator, including secondary electrons energy deposit 144 DtotalLetD[voxel] += DE + DEEletrons; // total dose Let Denominator, including secondary electrons energy deposit 145 totalLetT[voxel] += DX * Lsn; // total track Let Numerator 146 DtotalLetT[voxel] += DX; // total track Let Denominator 147 148 // store primary ions and secondary ions 149 size_t l; 150 for (l=0; l < ionLetStore.size(); l++) 151 { 152 // judge species of the current particle and store it 153 if (ionLetStore[l].PDGencoding == PDGencoding) 154 if ( ((trackID ==1) && (ionLetStore[l].isPrimary)) || ((trackID !=1) && (!ionLetStore[l].isPrimary))) 155 break; 156 } 157 158 if (l == ionLetStore.size()) // Just another type of ion/particle for our store... 159 { 160 // mass number 161 G4int A = particleDef -> GetAtomicMass(); 162 163 // particle name 164 G4String fullName = particleDef -> GetParticleName(); 165 G4String name = fullName.substr (0, fullName.find("[") ); // cut excitation energy [x.y] 166 167 ionLet ion = 168 { 169 (trackID == 1) ? true:false, // is it the primary particle? 170 PDGencoding, 171 fullName, 172 name, 173 Z, 174 A, 175 new G4double[nVoxels], // Let Dose Numerator 176 new G4double[nVoxels], // Let Dose Denominator 177 new G4double[nVoxels], // Let Track Numerator 178 new G4double[nVoxels], // Let Track Denominator 179 }; 180 181 // Initialize let and other parameters 182 for(G4int v=0; v < nVoxels; v++) 183 { 184 ion.letDN[v] = ion.letDD[v] = ion.letTN[v] = ion.letTD[v] = 0.; 185 } 186 187 188 ionLetStore.push_back(ion); 189 } 190 191 // calculate ions let and store them 192 ionLetStore[l].letDN[voxel] += (DE + DEEletrons)* Lsn; // ions dose Let Numerator, including secondary electrons energy deposit 193 ionLetStore[l].letDD[voxel] += DE + DEEletrons; // ions dose Let Denominator, including secondary electrons energy deposit 194 ionLetStore[l].letTN[voxel] += DX* Lsn; // ions track Let Numerator 195 ionLetStore[l].letTD[voxel] += DX; // ions track Let Denominator 196 197 } 198 199 200 201 202 void HadrontherapyLet::LetOutput() 203 { 204 for(G4int v=0; v < nVoxels; v++) 205 { 206 // compute total let 207 if (DtotalLetD[v]>0.) totalLetD[v] = totalLetD[v]/DtotalLetD[v]; 208 if (DtotalLetT[v]>0.) totalLetT[v] = totalLetT[v]/DtotalLetT[v]; 209 } 210 211 // Sort ions by A and then by Z ... 212 std::sort(ionLetStore.begin(), ionLetStore.end()); 213 214 215 // Compute Let Track and Let Dose for ions 216 217 for(G4int v=0; v < nVoxels; v++) 218 { 219 220 for (size_t ion=0; ion < ionLetStore.size(); ion++) 221 { 222 // compute ions let 223 if (ionLetStore[ion].letDD[v] >0.) ionLetStore[ion].letDN[v] = ionLetStore[ion].letDN[v] / ionLetStore[ion].letDD[v]; 224 if (ionLetStore[ion].letTD[v] >0.) ionLetStore[ion].letTN[v] = ionLetStore[ion].letTN[v] / ionLetStore[ion].letTD[v]; 225 } // end loop over ions 226 } 227 } // end loop over voxels 228 229 230 231 void HadrontherapyLet::StoreLetAscii() 232 { 233 #define width 15L 234 235 if(ionLetStore.size()) 236 { 237 ofs.open(filename, std::ios::out); 238 if (ofs.is_open()) 239 { 240 241 // Write the voxels index and total Lets and the list of particles/ions 242 ofs << "i" << '\t' << "j" << '\t' << "k"; 243 244 ofs << '\t' << "LDT"; 245 ofs << '\t' << "LTT"; 246 247 for (size_t l=0; l < ionLetStore.size(); l++) // write ions name 248 { 249 G4String a = (ionLetStore[l].isPrimary) ? "_1_D":"_D"; 250 ofs << '\t' << ionLetStore[l].name + a ; 251 G4String b = (ionLetStore[l].isPrimary) ? "_1_T":"_T"; 252 ofs << '\t' << ionLetStore[l].name + b ; 253 } 254 255 256 // Write data 257 258 G4AnalysisManager* LetFragmentTuple = G4AnalysisManager::Instance(); 259 260 LetFragmentTuple->SetVerboseLevel(1); 261 LetFragmentTuple->SetFirstHistoId(1); 262 LetFragmentTuple->SetFirstNtupleId(1); 263 LetFragmentTuple ->OpenFile("Let.csv"); 264 265 266 LetFragmentTuple ->CreateNtuple("coordinate", "Let"); 267 268 269 LetFragmentTuple ->CreateNtupleIColumn("i");//0 270 LetFragmentTuple ->CreateNtupleIColumn("j");//1 271 LetFragmentTuple ->CreateNtupleIColumn("k");//2 272 LetFragmentTuple ->CreateNtupleDColumn("TotalLetD");//3 273 LetFragmentTuple ->CreateNtupleDColumn("TotalLetT");//4 274 LetFragmentTuple ->CreateNtupleIColumn("A");//5 275 LetFragmentTuple ->CreateNtupleIColumn("Z");//6 276 LetFragmentTuple ->CreateNtupleDColumn("IonLETD");//7 277 LetFragmentTuple ->CreateNtupleDColumn("IonLETT");//8 278 LetFragmentTuple ->FinishNtuple(); 279 280 281 for(G4int i = 0; i < numberOfVoxelAlongX; i++) 282 for(G4int j = 0; j < numberOfVoxelAlongY; j++) 283 for(G4int k = 0; k < numberOfVoxelAlongZ; k++) 284 { 285 LetFragmentTuple->FillNtupleIColumn(1,0, i); 286 LetFragmentTuple->FillNtupleIColumn(1,1, j); 287 LetFragmentTuple->FillNtupleIColumn(1,2, k); 288 289 G4int v = matrix -> Index(i, j, k); 290 291 // write total Lets and voxels index 292 ofs << G4endl; 293 ofs << i << '\t' << j << '\t' << k ; 294 ofs << '\t' << totalLetD[v]/(keV/um); 295 ofs << '\t' << totalLetT[v]/(keV/um); 296 297 298 // write ions Lets 299 for (size_t l=0; l < ionLetStore.size(); l++) 300 { 301 302 // Write only not identically null data line 303 if(ionLetStore[l].letDN) 304 { 305 // write ions Lets 306 ofs << '\t' << ionLetStore[l].letDN[v]/(keV/um) ; 307 ofs << '\t' << ionLetStore[l].letTN[v]/(keV/um); 308 } 309 } 310 311 LetFragmentTuple->FillNtupleDColumn(1,3, totalLetD[v]/(keV/um)); 312 LetFragmentTuple->FillNtupleDColumn(1,4, totalLetT[v]/(keV/um)); 313 314 315 for (size_t ll=0; ll < ionLetStore.size(); ll++) 316 { 317 318 319 LetFragmentTuple->FillNtupleIColumn(1,5, ionLetStore[ll].A); 320 LetFragmentTuple->FillNtupleIColumn(1,6, ionLetStore[ll].Z); 321 322 323 LetFragmentTuple->FillNtupleDColumn(1,7, ionLetStore[ll].letDN[v]/(keV/um)); 324 LetFragmentTuple->FillNtupleDColumn(1,8, ionLetStore[ll].letTN[v]/(keV/um)); 325 LetFragmentTuple->AddNtupleRow(1); 326 } 327 } 328 ofs.close(); 329 330 LetFragmentTuple->Write(); 331 LetFragmentTuple->CloseFile(); 332 } 333 334 } 335 336 } 337 338