| 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 // This example is provided by the Geant4-DNA collaboration
27 // dnadamage3 example is derived from the chem6 example
28 // chem6 example authors: W. G. Shin and S. Incerti (CENBG, France)
29 //
30 // Any report or published results obtained using the Geant4-DNA software
31 // shall cite the following Geant4-DNA collaboration publication:
32 // J. Appl. Phys. 125 (2019) 104301
33 // Med. Phys. 45 (2018) e722-e739
34 // J. Comput. Phys. 274 (2014) 841-882
35 // Med. Phys. 37 (2010) 4692-4708
36 // Int. J. Model. Simul. Sci. Comput. 1 (2010) 157-178
37 // The Geant4-DNA web site is available at http://geant4-dna.org
38 //
39 // Authors: J. Naoki D. Kondo (UCSF, US)
40 // J. Ramos-Mendez and B. Faddegon (UCSF, US)
41 //
42 /// \file ScoreSpecies.cc
43 /// \brief Implementation of the ScoreSpecies class
44 ///
45 /// This is a primitive scorer class for molecular species.
46 /// The number of species is recorded for all times (predetermined or
47 /// user chosen). It also scores the energy deposition in order to compute the
48 /// radiochemical yields.
49
50 #include "ScoreSpecies.hh"
51
52 #include "G4AnalysisManager.hh"
53 #include "G4Event.hh"
54 #include "G4Scheduler.hh"
55 #include "G4TScoreNtupleWriter.hh"
56 #include "G4UImessenger.hh"
57 #include "G4UnitsTable.hh"
58
59 #include <G4EventManager.hh>
60 #include <G4MolecularConfiguration.hh>
61 #include <G4MoleculeCounter.hh>
62 #include <G4SystemOfUnits.hh>
63 #include <globals.hh>
64
65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
66
67 ScoreSpecies::ScoreSpecies(G4String name, G4int depth)
68 : G4VPrimitiveScorer(name, depth),
69 G4UImessenger(),
70 fOutputType("root"), // other options: "csv", "hdf5", "xml"
71 fHCID(-1)
72 {
73 fSpeciesdir = new G4UIdirectory("/scorer/species/");
74 fSpeciesdir->SetGuidance("ScoreSpecies commands");
75
76 fAddTimeToRecordcmd = new G4UIcmdWithADoubleAndUnit("/scorer/species/addTimeToRecord", this);
77
78 fTimeBincmd = new G4UIcmdWithAnInteger("/scorer/species/nOfTimeBins", this);
79
80 fOutputTypeUI = new G4UIcmdWithAString("/scorer/species/OutputFormat", this);
81 fOutputFileUI = new G4UIcmdWithAString("/scorer/species/OutputFile", this);
82
83 G4MoleculeCounter::Instance()->ResetCounter();
84 }
85
86 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
87
88 ScoreSpecies::~ScoreSpecies()
89 {
90 delete fSpeciesdir;
91 delete fAddTimeToRecordcmd;
92 delete fTimeBincmd;
93 delete fOutputTypeUI;
94 delete fOutputFileUI;
95 }
96
97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
98
99 void ScoreSpecies::SetNewValue(G4UIcommand* command, G4String newValue)
100 {
101 if (command == fAddTimeToRecordcmd) {
102 G4double cmdTime = fAddTimeToRecordcmd->GetNewDoubleValue(newValue);
103 AddTimeToRecord(cmdTime);
104 }
105 if (command == fTimeBincmd) {
106 ClearTimeToRecord();
107 G4int cmdBins = fTimeBincmd->GetNewIntValue(newValue);
108 G4double timeMin = 1 * ps;
109 G4double timeMax = G4Scheduler::Instance()->GetEndTime() - 1 * ps;
110 G4double timeLogMin = std::log10(timeMin);
111 G4double timeLogMax = std::log10(timeMax);
112 for (G4int i = 0; i < cmdBins; i++) {
113 AddTimeToRecord(std::pow(10, timeLogMin + i * (timeLogMax - timeLogMin) / (cmdBins - 1)));
114 }
115 }
116
117 if (command == fOutputTypeUI) {
118 fOutputType = newValue;
119 G4StrUtil::to_lower(fOutputType);
120 }
121
122 if (command == fOutputFileUI) fOutputFile = newValue;
123 }
124
125 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
126
127 G4bool ScoreSpecies::ProcessHits(G4Step* aStep, G4TouchableHistory*)
128 {
129 G4double edep = aStep->GetTotalEnergyDeposit();
130
131 if (edep == 0.) return FALSE;
132
133 edep *= aStep->GetPreStepPoint()->GetWeight();
134 G4int index = GetIndex(aStep);
135 fEvtMap->add(index, edep);
136 fEdep += edep;
137
138 return TRUE;
139 }
140
141 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
142
143 void ScoreSpecies::Initialize(G4HCofThisEvent* HCE)
144 {
145 fEvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName());
146
147 if (fHCID < 0) {
148 fHCID = GetCollectionID(0);
149 }
150
151 HCE->AddHitsCollection(fHCID, (G4VHitsCollection*)fEvtMap);
152 G4MoleculeCounter::Instance()->ResetCounter();
153 }
154
155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
156
157 void ScoreSpecies::EndOfEvent(G4HCofThisEvent*)
158 {
159 if (G4EventManager::GetEventManager()->GetConstCurrentEvent()->IsAborted()) {
160 fEdep = 0.;
161 G4MoleculeCounter::Instance()->ResetCounter();
162 return;
163 }
164
165 auto species = G4MoleculeCounter::Instance()->GetRecordedMolecules();
166
167 if (species.get() == 0 || species->size() == 0) {
168 G4cout << "No molecule recorded, energy deposited= " << G4BestUnit(fEdep, "Energy") << G4endl;
169 ++fNEvent;
170 fEdep = 0.;
171 G4MoleculeCounter::Instance()->ResetCounter();
172 return;
173 }
174
175 for (auto molecule : *species) {
176 for (auto time_mol : fTimeToRecord) {
177 double n_mol = G4MoleculeCounter::Instance()->GetNMoleculesAtTime(molecule, time_mol);
178
179 if (n_mol < 0) {
180 G4cerr << "N molecules not valid < 0 " << G4endl;
181 G4Exception("", "N<0", FatalException, "");
182 }
183
184 SpeciesInfo& molInfo = fSpeciesInfoPerTime[time_mol][molecule];
185 molInfo.fNumber += n_mol;
186 molInfo.fNumber2 += n_mol * n_mol;
187 double gValue = (n_mol / (fEdep / eV)) * 100.;
188 molInfo.fG += gValue;
189 molInfo.fG2 += gValue * gValue;
190 }
191 }
192
193 ++fNEvent;
194
195 fEdep = 0.;
196 G4MoleculeCounter::Instance()->ResetCounter();
197 }
198
199 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
200
201 void ScoreSpecies::AbsorbResultsFromWorkerScorer(G4VPrimitiveScorer* workerScorer)
202 {
203 ScoreSpecies* right = dynamic_cast<ScoreSpecies*>(workerScorer);
204
205 if (right == 0) {
206 return;
207 }
208 if (right == this) {
209 return;
210 }
211
212 SpeciesMap::iterator it_map1 = right->fSpeciesInfoPerTime.begin();
213 SpeciesMap::iterator end_map1 = right->fSpeciesInfoPerTime.end();
214
215 for (; it_map1 != end_map1; ++it_map1) {
216 InnerSpeciesMap& map2 = it_map1->second;
217 InnerSpeciesMap::iterator it_map2 = map2.begin();
218 InnerSpeciesMap::iterator end_map2 = map2.end();
219
220 for (; it_map2 != end_map2; ++it_map2) {
221 SpeciesInfo& molInfo = fSpeciesInfoPerTime[it_map1->first][it_map2->first];
222 molInfo.fNumber += it_map2->second.fNumber;
223 molInfo.fNumber2 += it_map2->second.fNumber2;
224 molInfo.fG += it_map2->second.fG;
225 molInfo.fG2 += it_map2->second.fG2;
226 }
227 }
228 right->fSpeciesInfoPerTime.clear();
229
230 fNEvent += right->fNEvent;
231 right->fNEvent = 0;
232 right->fEdep = 0.;
233 }
234
235 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
236
237 void ScoreSpecies::DrawAll() {}
238
239 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
240
241 void ScoreSpecies::PrintAll()
242 {
243 G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
244 G4cout << " PrimitiveScorer " << GetName() << G4endl;
245 G4cout << " Number of events " << fNEvent << G4endl;
246 G4cout << " Number of energy deposition recorded " << fEvtMap->entries() << G4endl;
247
248 for (auto itr : *fEvtMap->GetMap()) {
249 G4cout << " copy no.: " << itr.first << " energy deposit: " << *(itr.second) / GetUnitValue()
250 << " [" << GetUnit() << "]" << G4endl;
251 }
252 }
253
254 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
255
256 void ScoreSpecies::OutputAndClear()
257 {
258 if (G4Threading::IsWorkerThread()) return;
259
260 //---------------------------------------------------------------------------
261 // Save results
262
263 if (fOutputType != "ascii") {
264 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
265 analysisManager->SetDefaultFileType(fOutputType);
266
267 if (analysisManager) {
268 this->WriteWithAnalysisManager(analysisManager);
269 }
270 }
271 else
272 OutputToASCII();
273
274 fNEvent = 0;
275 fSpeciesInfoPerTime.clear();
276 }
277
278 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
279
280 void ScoreSpecies::WriteWithAnalysisManager(G4VAnalysisManager* analysisManager)
281 {
282 analysisManager->OpenFile(fOutputFile);
283 int fNtupleID = analysisManager->CreateNtuple("species", "species");
284 analysisManager->CreateNtupleIColumn(fNtupleID, "speciesID");
285 analysisManager->CreateNtupleIColumn(fNtupleID, "number");
286 analysisManager->CreateNtupleIColumn(fNtupleID, "nEvent");
287 analysisManager->CreateNtupleSColumn(fNtupleID, "speciesName");
288 analysisManager->CreateNtupleDColumn(fNtupleID, "time");
289 analysisManager->CreateNtupleDColumn(fNtupleID, "sumG");
290 analysisManager->CreateNtupleDColumn(fNtupleID, "sumG2");
291 analysisManager->FinishNtuple(fNtupleID);
292
293 for (auto it_map1 : fSpeciesInfoPerTime) {
294 InnerSpeciesMap& map2 = it_map1.second;
295
296 for (auto it_map2 : map2) {
297 double time = it_map1.first;
298 auto species = it_map2.first;
299 const G4String& name = species->GetName();
300 int molID = it_map2.first->GetMoleculeID();
301 int number = it_map2.second.fNumber;
302 double G = it_map2.second.fG;
303 double G2 = it_map2.second.fG2;
304
305 analysisManager->FillNtupleIColumn(fNtupleID, 0, molID);
306 analysisManager->FillNtupleIColumn(fNtupleID, 1, number);
307 analysisManager->FillNtupleIColumn(fNtupleID, 2, fNEvent);
308 analysisManager->FillNtupleSColumn(fNtupleID, 3, name);
309 analysisManager->FillNtupleDColumn(fNtupleID, 4, time);
310 analysisManager->FillNtupleDColumn(fNtupleID, 5, G);
311 analysisManager->FillNtupleDColumn(fNtupleID, 6, G2);
312 analysisManager->AddNtupleRow(fNtupleID);
313 }
314 }
315
316 analysisManager->Write();
317 analysisManager->CloseFile();
318 fRunID++;
319 }
320
321 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
322
323 void ScoreSpecies::OutputToASCII()
324 {
325 std::ofstream SpeciesOutput;
326 SpeciesOutput.open(fOutputFile + ".txt");
327
328 std::map<G4String, std::map<G4double, std::vector<G4double>>> mol;
329
330 for (auto it_map1 : fSpeciesInfoPerTime) {
331 InnerSpeciesMap& map2 = it_map1.second;
332 G4double time = it_map1.first / ps;
333
334 for (auto it_map2 : map2) {
335 G4double G = it_map2.second.fG;
336 G4double G2 = it_map2.second.fG2;
337 G4double Nb = it_map2.second.fNumber;
338 G4double Nb2 = it_map2.second.fNumber2;
339 G4double N = fNEvent;
340 if (N > 0) {
341 G /= N;
342 Nb /= N;
343 }
344
345 if (N == 1) {
346 G2 = 0.0;
347 Nb2 = 0.0;
348 }
349
350 else if (N > 0) {
351 G2 = std::sqrt(N / (N - 1) * (G2 / N - G * G));
352 Nb2 = std::sqrt(N / (N - 1) * (Nb2 / N - Nb * Nb));
353 }
354
355 mol[it_map2.first->GetName()][time].push_back(G);
356 mol[it_map2.first->GetName()][time].push_back(G2);
357 mol[it_map2.first->GetName()][time].push_back(Nb);
358 mol[it_map2.first->GetName()][time].push_back(Nb2);
359 }
360 }
361
362 SpeciesOutput << "# Species Scorer Output " << G4endl;
363 SpeciesOutput << "# " << std::setw(10) << "Time [ps]" << std::setw(15) << "Gx(/100 eV)"
364 << std::setw(15) << "RMS" << std::setw(15) << "Gx(#Mols)" << std::setw(20)
365 << "Molecule" << G4endl;
366
367 for (auto it1 : mol)
368 for (auto it2 : it1.second)
369 SpeciesOutput << std::setw(12) << it2.first << " " << std::setw(12) << it2.second[0]
370 << " " << std::setw(12) << it2.second[1] << " " << std::setw(12)
371 << it2.second[2] << " " << std::setw(17) << it1.first << G4endl;
372
373 SpeciesOutput.close();
374 }
375
376 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
377