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