| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 /// \file radiobiology/src/LETAccumulable.cc
28 /// \brief Implementation of the RadioBio::LET
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........oo
47
48 LETAccumulable::LETAccumulable() : VRadiobiolo
49
50 //....oooOO0OOooo........oooOO0OOooo........oo
51
52 void LETAccumulable::Merge(const G4VAccumulabl
53 {
54 if (GetVerboseLevel() > 0) {
55 G4cout << "LETAccumulable::Merge()" << G4e
56 }
57 const LETAccumulable& other = dynamic_cast<c
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.
70 G4int PDGencoding = rhsIonLetStore[l].GetP
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() ==
76 if (rhsIonLetStore[l].IsPrimary() == f
77 }
78 }
79
80 if (q == fIonLetStore.size()) // Just cop
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........oo
88
89 void LETAccumulable::Reset()
90 {
91 if (GetVerboseLevel() > 0) {
92 G4cout << "LETAccumulable::Reset()" << G4e
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........oo
107
108 // To accumulate given the hit
109 void LETAccumulable::Accumulate(Hit* hit)
110 {
111 if (GetVerboseLevel() > 1) {
112 G4cout << "LETAccumulable::Accumulate()" <
113 }
114
115 // No need to accumulate if no energy deposi
116 if (hit->GetDeltaE() == 0. || hit->GetTrackL
117
118 // Atomic number
119 G4int Z = hit->GetPartType()->GetAtomicNumbe
120 if (Z < 1) return; // Calculate only proton
121
122 G4int PDGencoding = hit->GetPartType()->GetP
123
124 if (PDGencoding == 22 || PDGencoding == 11)
125
126 // Get simple Particle data Group ID without
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()-
138
139 // Get mean kinetic energy
140 G4double ekinMean = hit->GetEkinMean();
141
142 // Get particle definition
143 const G4ParticleDefinition* particleDef = hi
144
145 // Get material
146 G4Material* mat = hit->GetPhysicalVolume()->
147
148 // ICRU stopping power calculation
149 G4EmCalculator emCal;
150 // Use the mean kinetic energy of ions in a
151 G4double Lsn = emCal.ComputeElectronicDEDX(e
152 // G4cout << "ERASE ME. Lsn = " << Lsn << G4
153 // G4cout << "ERASE ME. ekinMean = " << ekin
154
155 // Get deposited energy
156 G4double DE = hit->GetDeltaE();
157
158 // Get deposited energy considering secondar
159 G4double DEElectrons = hit->GetElectronEnerg
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 secon
169 fTotalLETD[voxel] += (DE + DEElectrons) * Ls
170 // Total dose LET Denominator, including sec
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 a
180 if (fIonLetStore[l].GetPDGencoding() == PD
181 if (((trackID == 1) && (fIonLetStore[l].
182 || ((trackID != 1) && (!fIonLetStore
183 break;
184 }
185
186 // Just another type of ion/particle for our
187 if (l == fIonLetStore.size()) {
188 // Mass number
189 G4int A = particleDef->GetAtomicMass();
190
191 // Particle name
192 G4String fullName = particleDef->GetPartic
193 G4String name = fullName.substr(0, fullNam
194
195 IonLet ion(trackID, PDGencoding, fullName,
196 VoxelizedSensitiveDetector::Get
197 fIonLetStore.push_back(ion);
198 }
199
200 // Calculate ions LET and store them
201 fIonLetStore[l].Update(voxel, DE, DEElectron
202 }
203
204 //....oooOO0OOooo........oooOO0OOooo........oo
205
206 G4int LETAccumulable::GetVerboseLevel() const
207 {
208 // return same level of LET class
209 return Manager::GetInstance()->GetQuantity("
210 }
211
212 //....oooOO0OOooo........oooOO0OOooo........oo
213
214 void LETAccumulable::Initialize()
215 {
216 if (GetVerboseLevel() > 0) {
217 G4cout << "LETAccumulable::Initialize(): "
218 }
219
220 G4int voxNumber = VoxelizedSensitiveDetector
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........oo
230
231 } // namespace RadioBio