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