| 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 // Any report or published results obtained using the Geant4-DNA software
28 // shall cite the following Geant4-DNA collaboration publication:
29 // Med. Phys. 37 (2010) 4692-4708
30 // J. Comput. Phys. 274 (2014) 841-882
31 // Phys. Med. Biol. 63(10) (2018) 105014-12pp
32 // The Geant4-DNA web site is available at http://geant4-dna.org
33 //
34 // ScoreSpecies.cc
35 //
36 #include "ScoreSpecies.hh"
37
38 #include "G4Event.hh"
39 #include "G4UnitsTable.hh"
40
41 #include <G4EventManager.hh>
42 #include <G4MolecularConfiguration.hh>
43 #include <G4MoleculeCounter.hh>
44 #include <G4SystemOfUnits.hh>
45 #include <globals.hh>
46 #include <iomanip>
47
48 /**
49 \file ScoreSpecies.cc
50 \class ScoreSpecies
51 This is a primitive scorer class for molecular species.
52 The number of species is recorded for all times (log spaced). It
53 also scores the energy deposition in order to compute the
54 radiochemical yields.
55 */
56
57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
58
59 ScoreSpecies::ScoreSpecies(G4String name, G4int depth)
60 : G4VPrimitiveScorer(name, depth), fEdep(0), fHCID(-1), fEvtMap(0)
61 {
62 fNEvent = 0;
63 G4double tMin = 1.0 * CLHEP::picosecond;
64 G4double tMax = 999999 * CLHEP::picosecond;
65 G4double tLogMin = std::log10(tMin);
66 G4double tLogMax = std::log10(tMax);
67 G4int tBins = 50;
68 for (int i = 0; i <= tBins; i++)
69 AddTimeToRecord(std::pow(10., tLogMin + i * (tLogMax - tLogMin) / tBins));
70
71 fEdep = 0;
72 }
73
74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
75
76 ScoreSpecies::~ScoreSpecies()
77 {
78 ;
79 }
80
81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
82
83 G4bool ScoreSpecies::ProcessHits(G4Step* aStep, G4TouchableHistory*)
84 {
85 G4double edep = aStep->GetTotalEnergyDeposit();
86
87 if (edep == 0.) return FALSE;
88
89 edep *= aStep->GetPreStepPoint()->GetWeight();
90 G4int index = GetIndex(aStep);
91 fEvtMap->add(index, edep);
92 fEdep += edep;
93
94 return TRUE;
95 }
96
97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
98
99 void ScoreSpecies::Initialize(G4HCofThisEvent* HCE)
100 {
101 fEvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName());
102
103 if (fHCID < 0) {
104 fHCID = GetCollectionID(0);
105 }
106
107 HCE->AddHitsCollection(fHCID, (G4VHitsCollection*)fEvtMap);
108 }
109
110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
111
112 void ScoreSpecies::EndOfEvent(G4HCofThisEvent*)
113 {
114 if (G4EventManager::GetEventManager()->GetConstCurrentEvent()->IsAborted()) {
115 fEdep = 0.;
116 G4MoleculeCounter::Instance()->ResetCounter();
117 return;
118 }
119
120 auto species = G4MoleculeCounter::Instance()->GetRecordedMolecules();
121
122 if (species.get() == 0 || species->size() == 0) {
123 G4cout << "No molecule recorded, energy deposited= " << G4BestUnit(fEdep, "Energy") << G4endl;
124 ++fNEvent;
125 fEdep = 0.;
126 G4MoleculeCounter::Instance()->ResetCounter();
127 return;
128 }
129
130 for (auto molecule : *species) {
131 for (auto time_mol : fTimeToRecord) {
132 double n_mol = G4MoleculeCounter::Instance()->GetNMoleculesAtTime(molecule, time_mol);
133
134 if (n_mol < 0) {
135 G4cerr << "N molecules not valid < 0 " << G4endl;
136 G4Exception("", "N<0", FatalException, "");
137 }
138
139 SpeciesInfo& molInfo = fSpeciesInfoPerTime[time_mol][molecule];
140 molInfo.fNumber += n_mol;
141 double gValue = (n_mol / (fEdep / eV)) * 100.;
142 molInfo.fG += gValue;
143 molInfo.fG2 += gValue * gValue;
144 }
145 }
146 ++fNEvent;
147
148 fEdep = 0.;
149 G4MoleculeCounter::Instance()->ResetCounter();
150 }
151
152 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
153
154 void ScoreSpecies::AbsorbResultsFromWorkerScorer(G4VPrimitiveScorer* workerScorer)
155 {
156 ScoreSpecies* right =
157 dynamic_cast<ScoreSpecies*>(dynamic_cast<G4VPrimitiveScorer*>(workerScorer));
158
159 if (right == 0) {
160 return;
161 }
162 if (right == this) {
163 return;
164 }
165
166 SpeciesMap::iterator it_map1 = right->fSpeciesInfoPerTime.begin();
167 SpeciesMap::iterator end_map1 = right->fSpeciesInfoPerTime.end();
168
169 for (; it_map1 != end_map1; ++it_map1) {
170 InnerSpeciesMap& map2 = it_map1->second;
171 InnerSpeciesMap::iterator it_map2 = map2.begin();
172 InnerSpeciesMap::iterator end_map2 = map2.end();
173
174 for (; it_map2 != end_map2; ++it_map2) {
175 SpeciesInfo& molInfo = fSpeciesInfoPerTime[it_map1->first][it_map2->first];
176 molInfo.fNumber += it_map2->second.fNumber;
177 molInfo.fG += it_map2->second.fG;
178 molInfo.fG2 += it_map2->second.fG2;
179 }
180 }
181
182 fNEvent += right->fNEvent;
183 right->fNEvent = 0;
184 right->fEdep = 0.;
185 }
186
187 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
188
189 void ScoreSpecies::DrawAll()
190 {
191 ;
192 }
193
194 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
195
196 void ScoreSpecies::PrintAll()
197 {
198 G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
199 G4cout << " PrimitiveScorer " << GetName() << G4endl;
200 G4cout << " Number of events " << fNEvent << G4endl;
201 G4cout << " Number of energy deposition recorded " << fEvtMap->entries() << G4endl;
202
203 for (auto itr : *fEvtMap->GetMap()) {
204 G4cout << " copy no.: " << itr.first << " energy deposit: " << *(itr.second) / GetUnitValue()
205 << " [" << GetUnit() << "]" << G4endl;
206 }
207 }
208
209 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
210
211 void ScoreSpecies::ASCII()
212 {
213 std::ofstream out("Species.txt");
214 if (!out) return;
215
216 out << "# Time [ps] G-value (/100 eV) RMS Molecule" << G4endl;
217
218 std::map<G4String, std::map<G4double, std::pair<G4double, G4double>>> mol;
219
220 for (auto it_map1 : fSpeciesInfoPerTime) {
221 InnerSpeciesMap& map2 = it_map1.second;
222 G4double time = it_map1.first / ps;
223 for (auto it_map2 : map2) {
224 G4double G = it_map2.second.fG;
225 G4double G2 = it_map2.second.fG2;
226 G4double N = fNEvent;
227 G /= N;
228 G2 = std::sqrt(N / (N - 1) * (G2 / N - G * G));
229 mol[it_map2.first->GetName()][time] = std::make_pair(G, G2);
230 }
231 }
232
233 for (auto it1 : mol)
234 for (auto it2 : it1.second)
235 out << std::setw(12) << it2.first << std::setw(12) << it2.second.first << std::setw(12)
236 << it2.second.second << std::setw(12) << std::setw(12) << it1.first << G4endl;
237
238 out.close();
239 }
240
241 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
242
243 void ScoreSpecies::OutputAndClear()
244 {
245 if (G4Threading::IsWorkerThread()) return;
246
247 //----------------------------------------------------------------------------
248 // Save results in ASCII format
249 ASCII();
250
251 fNEvent = 0;
252 fSpeciesInfoPerTime.clear();
253 }
254