| 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 //
27 /// \file radiobiology/src/LETAccumulable.cc
28 /// \brief Implementation of the RadioBio::LETAccumulable class
29
30 #include "LETAccumulable.hh"
31
32 #include "G4EmCalculator.hh"
33 #include "G4LogicalVolume.hh"
34 #include "G4ParticleDefinition.hh"
35
36 #include "Hit.hh"
37 #include "Manager.hh"
38 #include "VoxelizedSensitiveDetector.hh"
39 #include <G4SystemOfUnits.hh>
40
41 #include <tuple>
42
43 namespace RadioBio
44 {
45
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
47
48 LETAccumulable::LETAccumulable() : VRadiobiologicalAccumulable("LET") {}
49
50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
51
52 void LETAccumulable::Merge(const G4VAccumulable& rhs)
53 {
54 if (GetVerboseLevel() > 0) {
55 G4cout << "LETAccumulable::Merge()" << G4endl;
56 }
57 const LETAccumulable& other = dynamic_cast<const LETAccumulable&>(rhs);
58
59 // Merges standard counters
60 fTotalLETT += other.GetTotalLETT();
61 fDTotalLETT += other.GetDTotalLETT();
62 fDTotalLETD += other.GetDTotalLETD();
63 fTotalLETD += other.GetTotalLETD();
64
65 // Merges ion counters
66 auto rhsIonLetStore = other.GetIonLetStore();
67
68 // Loop over rhs ions
69 for (unsigned int l = 0; l < rhsIonLetStore.size(); l++) {
70 G4int PDGencoding = rhsIonLetStore[l].GetPDGencoding();
71 size_t q;
72 // Loop over lhs ions to find the one
73 for (q = 0; q < fIonLetStore.size(); q++) {
74 // Check if primary status is the same
75 if (fIonLetStore[q].GetPDGencoding() == PDGencoding) {
76 if (rhsIonLetStore[l].IsPrimary() == fIonLetStore[q].IsPrimary()) break;
77 }
78 }
79
80 if (q == fIonLetStore.size()) // Just copy the IonLet element in lhs store
81 fIonLetStore.push_back(rhsIonLetStore[l]);
82 else // Merges rhs data with lhs ones
83 fIonLetStore[q].Merge(&(rhsIonLetStore[l]));
84 }
85 }
86
87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
88
89 void LETAccumulable::Reset()
90 {
91 if (GetVerboseLevel() > 0) {
92 G4cout << "LETAccumulable::Reset()" << G4endl;
93 }
94 if (fInitialized) {
95 fTotalLETT = 0.0;
96 fDTotalLETT = 0.0;
97 fDTotalLETD = 0.0;
98 fTotalLETD = 0.0;
99 fIonLetStore.clear();
100 }
101 else {
102 Initialize();
103 }
104 }
105
106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
107
108 // To accumulate given the hit
109 void LETAccumulable::Accumulate(Hit* hit)
110 {
111 if (GetVerboseLevel() > 1) {
112 G4cout << "LETAccumulable::Accumulate()" << G4endl;
113 }
114
115 // No need to accumulate if no energy deposit or track length is zero.
116 if (hit->GetDeltaE() == 0. || hit->GetTrackLength() == 0.) return;
117
118 // Atomic number
119 G4int Z = hit->GetPartType()->GetAtomicNumber();
120 if (Z < 1) return; // Calculate only protons and ions
121
122 G4int PDGencoding = hit->GetPartType()->GetPDGEncoding();
123
124 if (PDGencoding == 22 || PDGencoding == 11) return; // do not accumulate for electrons or gammas
125
126 // Get simple Particle data Group ID without excitation level
127 PDGencoding -= PDGencoding % 10;
128
129 // Get voxel number
130 G4int xIndex = hit->GetXindex();
131 G4int yIndex = hit->GetYindex();
132 ;
133 G4int zIndex = hit->GetZindex();
134 ;
135
136 G4int voxel =
137 VoxelizedSensitiveDetector::GetInstance()->GetThisVoxelNumber(xIndex, yIndex, zIndex);
138
139 // Get mean kinetic energy
140 G4double ekinMean = hit->GetEkinMean();
141
142 // Get particle definition
143 const G4ParticleDefinition* particleDef = hit->GetPartType();
144
145 // Get material
146 G4Material* mat = hit->GetPhysicalVolume()->GetLogicalVolume()->GetMaterial();
147
148 // ICRU stopping power calculation
149 G4EmCalculator emCal;
150 // Use the mean kinetic energy of ions in a step to calculate ICRU stopping power
151 G4double Lsn = emCal.ComputeElectronicDEDX(ekinMean, particleDef, mat);
152 // G4cout << "ERASE ME. Lsn = " << Lsn << G4endl;
153 // G4cout << "ERASE ME. ekinMean = " << ekinMean << G4endl;
154
155 // Get deposited energy
156 G4double DE = hit->GetDeltaE();
157
158 // Get deposited energy considering secondary electrons
159 G4double DEElectrons = hit->GetElectronEnergy();
160
161 // Get track length
162 G4double DX = hit->GetTrackLength();
163
164 // Get track ID
165 G4int trackID = hit->GetTrackID();
166
167 // Total LET calculation...
168 // Total dose LET Numerator, including secondary electrons energy deposit
169 fTotalLETD[voxel] += (DE + DEElectrons) * Lsn;
170 // Total dose LET Denominator, including secondary electrons energy deposit
171 fDTotalLETD[voxel] += DE + DEElectrons;
172 // Total track LET Numerator
173 fTotalLETT[voxel] += DX * Lsn;
174 // Total track LET Denominator
175 fDTotalLETT[voxel] += DX;
176
177 size_t l;
178 for (l = 0; l < fIonLetStore.size(); l++) {
179 // Judge species of the current particle and store it
180 if (fIonLetStore[l].GetPDGencoding() == PDGencoding)
181 if (((trackID == 1) && (fIonLetStore[l].IsPrimary()))
182 || ((trackID != 1) && (!fIonLetStore[l].IsPrimary())))
183 break;
184 }
185
186 // Just another type of ion/particle for our store...
187 if (l == fIonLetStore.size()) {
188 // Mass number
189 G4int A = particleDef->GetAtomicMass();
190
191 // Particle name
192 G4String fullName = particleDef->GetParticleName();
193 G4String name = fullName.substr(0, fullName.find("[")); // Cut excitation energy [x.y]
194
195 IonLet ion(trackID, PDGencoding, fullName, name, Z, A,
196 VoxelizedSensitiveDetector::GetInstance()->GetTotalVoxelNumber());
197 fIonLetStore.push_back(ion);
198 }
199
200 // Calculate ions LET and store them
201 fIonLetStore[l].Update(voxel, DE, DEElectrons, Lsn, DX);
202 }
203
204 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
205
206 G4int LETAccumulable::GetVerboseLevel() const
207 {
208 // return same level of LET class
209 return Manager::GetInstance()->GetQuantity("LET")->GetVerboseLevel();
210 }
211
212 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
213
214 void LETAccumulable::Initialize()
215 {
216 if (GetVerboseLevel() > 0) {
217 G4cout << "LETAccumulable::Initialize(): " << G4endl;
218 }
219
220 G4int voxNumber = VoxelizedSensitiveDetector::GetInstance()->GetTotalVoxelNumber();
221
222 fTotalLETT = array_type(0.0, voxNumber);
223 fTotalLETD = array_type(0.0, voxNumber);
224 fDTotalLETT = array_type(0.0, voxNumber);
225 fDTotalLETD = array_type(0.0, voxNumber);
226 fInitialized = true;
227 }
228
229 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
230
231 } // namespace RadioBio