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