| 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/LET.cc
28 /// \brief Implementation of the RadioBio::LET class
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........oooOO0OOooo........oooOO0OOooo......
47
48 LET::LET() : VRadiobiologicalQuantity()
49 {
50 fPath = "RadioBio_LET.out"; // Default output filename
51
52 fMessenger = new LETMessenger(this);
53
54 Initialize();
55 }
56
57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
58
59 LET::~LET()
60 {
61 delete fMessenger;
62 }
63
64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
65
66 void LET::Initialize()
67 {
68 G4int VoxelNumber = VoxelizedSensitiveDetector::GetInstance()->GetTotalVoxelNumber();
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........oooOO0OOooo........oooOO0OOooo......
83
84 void LET::Compute()
85 {
86 // Skip LET computation if calculation not enabled.
87 if (!fCalculationEnabled) {
88 if (fVerboseLevel > 0) {
89 G4cout << "LET::Compute() called but skipped as calculation not enabled" << G4endl;
90 }
91 return;
92 }
93 if (fVerboseLevel > 0) G4cout << "LET::Compute()" << G4endl;
94
95 if (VoxelizedSensitiveDetector::GetInstance() == nullptr)
96 G4Exception("LET::ComputeLET", "PointerNotAvailable", FatalException,
97 "Computing LET without voxelized geometry pointer!");
98
99 G4int VoxelNumber = VoxelizedSensitiveDetector::GetInstance()->GetTotalVoxelNumber();
100
101 for (G4int v = 0; v < VoxelNumber; v++) {
102 if (fVerboseLevel > 1) G4cout << "COMPUTING LET of voxel " << v << G4endl;
103
104 // Compute total LET
105 if (fDTotalLETD[v] > 0.) fTotalLETD[v] = fNTotalLETD[v] / fDTotalLETD[v];
106 // G4cout << "ERASE ME. DEBUG: voxel = " << v << " fNTotalLETD[v] = " << fNTotalLETD[v]
107 // << " fDTotalLETD[v] = " << fDTotalLETD[v] << G4endl;
108 if (fDTotalLETT[v] > 0.) fTotalLETT[v] = fNTotalLETT[v] / fDTotalLETT[v];
109 }
110
111 // Sort ions by A and then by Z ...
112 std::sort(fIonLetStore.begin(), fIonLetStore.end());
113
114 // Compute LET Track and LET Dose for ions
115 for (size_t ion = 0; ion < fIonLetStore.size(); ion++) {
116 fIonLetStore[ion].Calculate();
117 }
118
119 fCalculated = true;
120 }
121
122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
123
124 // save LET
125 void LET::Store()
126 {
127 // Skip LET storing if calculation not enabled.
128 if (!fCalculationEnabled) {
129 if (fVerboseLevel > 0) {
130 G4cout << "LET::Store() called but skipped as calculation not enabled" << G4endl;
131 }
132 return;
133 }
134
135 if (fSaved == true)
136 G4Exception("LET::StoreLET", "NtuplesAlreadySaved", FatalException,
137 "Overwriting LET file. For multiple runs, change filename.");
138
139 Compute();
140 #define width 15L
141 if (fVerboseLevel > 0) G4cout << "LET::StoreLET" << G4endl;
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 and the list of particles/ions
149 ofs << std::setprecision(6) << std::left << "i\tj\tk\t";
150 ofs << std::setw(width) << "LDT";
151 ofs << std::setw(width) << "LTT";
152
153 for (size_t l = 0; l < fIonLetStore.size(); l++) // Write ions name
154 {
155 G4String a = (fIonLetStore[l].IsPrimary()) ? "_1_D" : "_D";
156 ofs << std::setw(width) << fIonLetStore[l].GetName() + a;
157 G4String b = (fIonLetStore[l].IsPrimary()) ? "_1_T" : "_T";
158 ofs << std::setw(width) << fIonLetStore[l].GetName() + b;
159 }
160 ofs << G4endl;
161
162 // Write data
163 G4AnalysisManager* LETFragmentTuple = G4AnalysisManager::Instance();
164
165 LETFragmentTuple->SetVerboseLevel(1);
166 LETFragmentTuple->SetFirstHistoId(1);
167 LETFragmentTuple->SetFirstNtupleId(1);
168
169 LETFragmentTuple->SetDefaultFileType("xml");
170 LETFragmentTuple->OpenFile("LET");
171
172 LETFragmentTuple->CreateNtuple("coordinate", "LET");
173
174 LETFragmentTuple->CreateNtupleIColumn("i"); // 0
175 LETFragmentTuple->CreateNtupleIColumn("j"); // 1
176 LETFragmentTuple->CreateNtupleIColumn("k"); // 2
177 LETFragmentTuple->CreateNtupleDColumn("TotalLETD"); // 3
178 LETFragmentTuple->CreateNtupleDColumn("TotalLETT"); // 4
179 LETFragmentTuple->CreateNtupleIColumn("A"); // 5
180 LETFragmentTuple->CreateNtupleIColumn("Z"); // 6
181 LETFragmentTuple->CreateNtupleDColumn("IonLetD"); // 7
182 LETFragmentTuple->CreateNtupleDColumn("IonLetT"); // 8
183 LETFragmentTuple->FinishNtuple();
184
185 auto voxSensDet = VoxelizedSensitiveDetector::GetInstance();
186
187 for (G4int i = 0; i < voxSensDet->GetVoxelNumberAlongX(); i++)
188 for (G4int j = 0; j < voxSensDet->GetVoxelNumberAlongY(); j++)
189 for (G4int k = 0; k < voxSensDet->GetVoxelNumberAlongZ(); k++) {
190 LETFragmentTuple->FillNtupleIColumn(1, 0, i);
191 LETFragmentTuple->FillNtupleIColumn(1, 1, j);
192 LETFragmentTuple->FillNtupleIColumn(1, 2, k);
193
194 G4int v = voxSensDet->GetThisVoxelNumber(i, j, k);
195
196 // Write total LETs and voxels index
197 ofs << G4endl;
198 ofs << i << '\t' << j << '\t' << k << '\t';
199 ofs << std::setw(width) << fTotalLETD[v] / (keV / um);
200 ofs << std::setw(width) << fTotalLETT[v] / (keV / um);
201
202 // Write ions LETs
203 for (size_t l = 0; l < fIonLetStore.size(); l++) {
204 // Write ions LETs
205 ofs << std::setw(width) << fIonLetStore[l].GetLETD()[v] / (keV / um);
206 ofs << std::setw(width) << fIonLetStore[l].GetLETT()[v] / (keV / um);
207 }
208
209 LETFragmentTuple->FillNtupleDColumn(1, 3, fTotalLETD[v] / (keV / um));
210 LETFragmentTuple->FillNtupleDColumn(1, 4, fTotalLETT[v] / (keV / um));
211
212 for (size_t ll = 0; ll < fIonLetStore.size(); ll++) {
213 LETFragmentTuple->FillNtupleIColumn(1, 5, fIonLetStore[ll].GetA());
214 LETFragmentTuple->FillNtupleIColumn(1, 6, fIonLetStore[ll].GetZ());
215
216 LETFragmentTuple->FillNtupleDColumn(1, 7, fIonLetStore[ll].GetLETD()[v] / (keV / um));
217 LETFragmentTuple->FillNtupleDColumn(1, 8, fIonLetStore[ll].GetLETT()[v] / (keV / um));
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........oooOO0OOooo........oooOO0OOooo......
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........oooOO0OOooo........oooOO0OOooo......
251
252 // Add data taken from accumulables
253 void LET::AddFromAccumulable(G4VAccumulable* GenAcc)
254 {
255 LETAccumulable* acc = (LETAccumulable*)GenAcc;
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->GetIonLetStore().size(); l++) {
265 G4int PDGencoding = acc->GetIonLetStore()[l].GetPDGencoding();
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() == PDGencoding) {
271 if (acc->GetIonLetStore()[l].IsPrimary() == fIonLetStore[q].IsPrimary()) break;
272 }
273 }
274 // If ion is missing, copy it
275 if (q == fIonLetStore.size())
276 fIonLetStore.push_back(acc->GetIonLetStore()[l]);
277 else // Merge rhs data with lhs ones
278 fIonLetStore[q].Merge(&(acc->GetIonLetStore()[l]));
279 }
280 fCalculated = false;
281 }
282
283 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
284
285 void LET::PrintParameters()
286 {
287 G4cout << "*******************************************" << G4endl
288 << "******* Parameters of the class LET *******" << G4endl
289 << "*******************************************" << G4endl;
290 PrintVirtualParameters();
291 }
292
293 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
294
295 } // namespace RadioBio