| 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/LET.cc
28 /// \brief Implementation of the RadioBio::LET
29
30 #include "LET.hh"
31
32 #include "G4EmCalculator.hh"
33 #include "G4RunManager.hh"
34 #include "G4SystemOfUnits.hh"
35
36 #include "DetectorConstruction.hh"
37 #include "Hit.hh"
38 #include "LETAccumulable.hh"
39 #include "LETMessenger.hh"
40 #include "VoxelizedSensitiveDetector.hh"
41
42 #include <cmath>
43
44 namespace RadioBio
45 {
46 //....oooOO0OOooo........oooOO0OOooo........oo
47
48 LET::LET() : VRadiobiologicalQuantity()
49 {
50 fPath = "RadioBio_LET.out"; // Default outp
51
52 fMessenger = new LETMessenger(this);
53
54 Initialize();
55 }
56
57 //....oooOO0OOooo........oooOO0OOooo........oo
58
59 LET::~LET()
60 {
61 delete fMessenger;
62 }
63
64 //....oooOO0OOooo........oooOO0OOooo........oo
65
66 void LET::Initialize()
67 {
68 G4int VoxelNumber = VoxelizedSensitiveDetect
69
70 // Instances for Total LET
71 fNTotalLETT.resize(VoxelNumber);
72 fNTotalLETD.resize(VoxelNumber);
73 fDTotalLETT.resize(VoxelNumber);
74 fDTotalLETD.resize(VoxelNumber);
75
76 fTotalLETD.resize(VoxelNumber);
77 fTotalLETT.resize(VoxelNumber);
78
79 fInitialized = true;
80 }
81
82 //....oooOO0OOooo........oooOO0OOooo........oo
83
84 void LET::Compute()
85 {
86 // Skip LET computation if calculation not e
87 if (!fCalculationEnabled) {
88 if (fVerboseLevel > 0) {
89 G4cout << "LET::Compute() called but ski
90 }
91 return;
92 }
93 if (fVerboseLevel > 0) G4cout << "LET::Compu
94
95 if (VoxelizedSensitiveDetector::GetInstance(
96 G4Exception("LET::ComputeLET", "PointerNot
97 "Computing LET without voxeliz
98
99 G4int VoxelNumber = VoxelizedSensitiveDetect
100
101 for (G4int v = 0; v < VoxelNumber; v++) {
102 if (fVerboseLevel > 1) G4cout << "COMPUTIN
103
104 // Compute total LET
105 if (fDTotalLETD[v] > 0.) fTotalLETD[v] = f
106 // G4cout << "ERASE ME. DEBUG: voxel = " <
107 // << " fDTotalLETD[v] = " << fDTot
108 if (fDTotalLETT[v] > 0.) fTotalLETT[v] = f
109 }
110
111 // Sort ions by A and then by Z ...
112 std::sort(fIonLetStore.begin(), fIonLetStore
113
114 // Compute LET Track and LET Dose for ions
115 for (size_t ion = 0; ion < fIonLetStore.size
116 fIonLetStore[ion].Calculate();
117 }
118
119 fCalculated = true;
120 }
121
122 //....oooOO0OOooo........oooOO0OOooo........oo
123
124 // save LET
125 void LET::Store()
126 {
127 // Skip LET storing if calculation not enabl
128 if (!fCalculationEnabled) {
129 if (fVerboseLevel > 0) {
130 G4cout << "LET::Store() called but skipp
131 }
132 return;
133 }
134
135 if (fSaved == true)
136 G4Exception("LET::StoreLET", "NtuplesAlrea
137 "Overwriting LET file. For mul
138
139 Compute();
140 #define width 15L
141 if (fVerboseLevel > 0) G4cout << "LET::Store
142
143 // If there is at least one ion
144 if (fIonLetStore.size()) {
145 std::ofstream ofs(fPath);
146 // ofs.open(fPath, std::ios::out);
147 if (ofs.is_open()) {
148 // Write the voxels index and total LETs
149 ofs << std::setprecision(6) << std::left
150 ofs << std::setw(width) << "LDT";
151 ofs << std::setw(width) << "LTT";
152
153 for (size_t l = 0; l < fIonLetStore.size
154 {
155 G4String a = (fIonLetStore[l].IsPrimar
156 ofs << std::setw(width) << fIonLetStor
157 G4String b = (fIonLetStore[l].IsPrimar
158 ofs << std::setw(width) << fIonLetStor
159 }
160 ofs << G4endl;
161
162 // Write data
163 G4AnalysisManager* LETFragmentTuple = G4
164
165 LETFragmentTuple->SetVerboseLevel(1);
166 LETFragmentTuple->SetFirstHistoId(1);
167 LETFragmentTuple->SetFirstNtupleId(1);
168
169 LETFragmentTuple->SetDefaultFileType("xm
170 LETFragmentTuple->OpenFile("LET");
171
172 LETFragmentTuple->CreateNtuple("coordina
173
174 LETFragmentTuple->CreateNtupleIColumn("i
175 LETFragmentTuple->CreateNtupleIColumn("j
176 LETFragmentTuple->CreateNtupleIColumn("k
177 LETFragmentTuple->CreateNtupleDColumn("T
178 LETFragmentTuple->CreateNtupleDColumn("T
179 LETFragmentTuple->CreateNtupleIColumn("A
180 LETFragmentTuple->CreateNtupleIColumn("Z
181 LETFragmentTuple->CreateNtupleDColumn("I
182 LETFragmentTuple->CreateNtupleDColumn("I
183 LETFragmentTuple->FinishNtuple();
184
185 auto voxSensDet = VoxelizedSensitiveDete
186
187 for (G4int i = 0; i < voxSensDet->GetVox
188 for (G4int j = 0; j < voxSensDet->GetV
189 for (G4int k = 0; k < voxSensDet->Ge
190 LETFragmentTuple->FillNtupleIColum
191 LETFragmentTuple->FillNtupleIColum
192 LETFragmentTuple->FillNtupleIColum
193
194 G4int v = voxSensDet->GetThisVoxel
195
196 // Write total LETs and voxels ind
197 ofs << G4endl;
198 ofs << i << '\t' << j << '\t' << k
199 ofs << std::setw(width) << fTotalL
200 ofs << std::setw(width) << fTotalL
201
202 // Write ions LETs
203 for (size_t l = 0; l < fIonLetStor
204 // Write ions LETs
205 ofs << std::setw(width) << fIonL
206 ofs << std::setw(width) << fIonL
207 }
208
209 LETFragmentTuple->FillNtupleDColum
210 LETFragmentTuple->FillNtupleDColum
211
212 for (size_t ll = 0; ll < fIonLetSt
213 LETFragmentTuple->FillNtupleICol
214 LETFragmentTuple->FillNtupleICol
215
216 LETFragmentTuple->FillNtupleDCol
217 LETFragmentTuple->FillNtupleDCol
218 LETFragmentTuple->AddNtupleRow(1
219 }
220 }
221 ofs.close();
222
223 // LETFragmentTuple->Write();
224 LETFragmentTuple->CloseFile();
225 }
226 }
227
228 fSaved = true;
229 }
230
231 //....oooOO0OOooo........oooOO0OOooo........oo
232
233 void LET::Reset()
234 {
235 if (fVerboseLevel > 1) {
236 G4cout << "LET::Reset(): ";
237 }
238 fNTotalLETT = 0.0;
239 fNTotalLETD = 0.0;
240 fDTotalLETT = 0.0;
241 fDTotalLETD = 0.0;
242
243 fTotalLETD = 0.0;
244 fTotalLETT = 0.0;
245
246 fIonLetStore.clear();
247 fCalculated = false;
248 }
249
250 //....oooOO0OOooo........oooOO0OOooo........oo
251
252 // Add data taken from accumulables
253 void LET::AddFromAccumulable(G4VAccumulable* G
254 {
255 LETAccumulable* acc = (LETAccumulable*)GenAc
256
257 AddNTotalLETT(acc->GetTotalLETT());
258 AddDTotalLETT(acc->GetDTotalLETT());
259 AddNTotalLETD(acc->GetTotalLETD());
260 AddDTotalLETD(acc->GetDTotalLETD());
261
262 // Merges ion counters
263 // Loop over rhs ions
264 for (unsigned int l = 0; l < acc->GetIonLetS
265 G4int PDGencoding = acc->GetIonLetStore()[
266 size_t q;
267 // Loop over lhs ions to find the one
268 for (q = 0; q < fIonLetStore.size(); q++)
269 // If he found it, sums values
270 if (fIonLetStore[q].GetPDGencoding() ==
271 if (acc->GetIonLetStore()[l].IsPrimary
272 }
273 }
274 // If ion is missing, copy it
275 if (q == fIonLetStore.size())
276 fIonLetStore.push_back(acc->GetIonLetSto
277 else // Merge rhs data with lhs ones
278 fIonLetStore[q].Merge(&(acc->GetIonLetSt
279 }
280 fCalculated = false;
281 }
282
283 //....oooOO0OOooo........oooOO0OOooo........oo
284
285 void LET::PrintParameters()
286 {
287 G4cout << "*********************************
288 << "******* Parameters of the class L
289 << "*********************************
290 PrintVirtualParameters();
291 }
292
293 //....oooOO0OOooo........oooOO0OOooo........oo
294
295 } // namespace RadioBio