| 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 // The Geant4-DNA web site is available at http://geant4-dna.org
32 //
33 // ScoreSpecies.cc
34 //
35 #include "ScoreSpecies.hh"
36
37 #include "G4Event.hh"
38 #include "G4UnitsTable.hh"
39
40 #include <G4AnalysisManager.hh>
41 #include <G4EventManager.hh>
42 #include <G4MolecularConfiguration.hh>
43 #include <G4MoleculeCounter.hh>
44 #include <G4SystemOfUnits.hh>
45 #include <globals.hh>
46
47 /**
48 \file ScoreSpecies.cc
49 \class ScoreSpecies
50 This is a primitive scorer class for molecular species.
51 The number of species is recorded for all times (predetermined or
52 user chosen). It also scores the energy deposition in order to compute the
53 radiochemical yields.
54 */
55
56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
57
58 ScoreSpecies::ScoreSpecies(G4String name, G4int depth)
59 : G4VPrimitiveScorer(name, depth),
60 fEdep(0),
61 fOutputType("root"), // other options: "csv", "hdf5", "xml"
62 fHCID(-1),
63 fEvtMap(0)
64 {
65 fNEvent = 0;
66 AddTimeToRecord(1 * CLHEP::picosecond);
67 AddTimeToRecord(10 * CLHEP::picosecond);
68 AddTimeToRecord(100 * CLHEP::picosecond);
69 AddTimeToRecord(1000 * CLHEP::picosecond);
70 AddTimeToRecord(10000 * CLHEP::picosecond);
71 AddTimeToRecord(100000 * CLHEP::picosecond);
72 AddTimeToRecord(999999 * CLHEP::picosecond);
73 fEdep = 0;
74 }
75
76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
77
78 ScoreSpecies::~ScoreSpecies()
79 {
80 ;
81 }
82
83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
84
85 G4bool ScoreSpecies::ProcessHits(G4Step* aStep, G4TouchableHistory*)
86 {
87 G4double edep = aStep->GetTotalEnergyDeposit();
88
89 if (edep == 0.) return FALSE;
90
91 edep *= aStep->GetPreStepPoint()->GetWeight(); // (Particle Weight)
92 G4int index = GetIndex(aStep);
93 fEvtMap->add(index, edep);
94 fEdep += edep;
95
96 return TRUE;
97 }
98
99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
100
101 void ScoreSpecies::Initialize(G4HCofThisEvent* HCE)
102 {
103 fEvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName());
104
105 if (fHCID < 0) {
106 fHCID = GetCollectionID(0);
107 }
108
109 HCE->AddHitsCollection(fHCID, (G4VHitsCollection*)fEvtMap);
110 }
111
112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
113
114 void ScoreSpecies::EndOfEvent(G4HCofThisEvent*)
115 {
116 if (G4EventManager::GetEventManager()->GetConstCurrentEvent()->IsAborted()) {
117 fEdep = 0.;
118 G4MoleculeCounter::Instance()->ResetCounter();
119 return;
120 }
121
122 auto species = G4MoleculeCounter::Instance()->GetRecordedMolecules();
123
124 if (species.get() == 0 || species->size() == 0) {
125 G4cout << "No molecule recorded, energy deposited= " << G4BestUnit(fEdep, "Energy") << G4endl;
126 ++fNEvent;
127 fEdep = 0.;
128 G4MoleculeCounter::Instance()->ResetCounter();
129 return;
130 }
131
132 // G4cout << "ScoreSpecies::EndOfEvent"<<G4endl;
133 // G4cout << "End of event, deposited energy: "
134 // << G4BestUnit(fEdep, "Energy") << G4endl;
135
136 #ifdef _ScoreSpecies_FOR_ALL_EVENTS
137 int eventID = G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetEventID();
138 #endif
139
140 for (auto molecule : *species) {
141 for (auto time_mol : fTimeToRecord) {
142 double n_mol = G4MoleculeCounter::Instance()->GetNMoleculesAtTime(molecule, time_mol);
143
144 if (n_mol < 0) {
145 G4cerr << "N molecules not valid < 0 " << G4endl;
146 G4Exception("", "N<0", FatalException, "");
147 }
148
149 SpeciesInfo& molInfo = fSpeciesInfoPerTime[time_mol][molecule];
150 molInfo.fNumber += n_mol;
151 double gValue = (n_mol / (fEdep / eV)) * 100.;
152 molInfo.fG += gValue;
153 molInfo.fG2 += gValue * gValue;
154
155 #ifdef _ScoreSpecies_FOR_ALL_EVENTS
156 SpeciesInfoSOA& molInfoPerEvent = fSpeciesInfoPerEvent[time_mol][molecule];
157 molInfoPerEvent.fNumber.push_back(n_mol);
158 molInfoPerEvent.fG.push_back(gValue);
159 molInfoPerEvent.fG2.push_back(gValue * gValue);
160 molInfoPerEvent.fEventID.push_back(eventID);
161 #endif
162 // G4cout << "In Save molucule: fNumber " << molInfo.fNumber
163 // << " fG " << molInfo.fG
164 // << " fEdep " << fEdep/eV
165 // << G4endl;
166 }
167 }
168
169 ++fNEvent;
170
171 // G4cout << "End of event " << fNEvent
172 // << ", energy deposited=" << G4BestUnit(fEdep, "Energy") << G4endl;
173
174 fEdep = 0.;
175 G4MoleculeCounter::Instance()->ResetCounter();
176 }
177
178 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
179
180 void ScoreSpecies::AbsorbResultsFromWorkerScorer(G4VPrimitiveScorer* workerScorer)
181 {
182 ScoreSpecies* right =
183 dynamic_cast<ScoreSpecies*>(dynamic_cast<G4VPrimitiveScorer*>(workerScorer));
184
185 if (right == 0) {
186 return;
187 }
188 if (right == this) {
189 return;
190 }
191
192 // G4cout<<"ScoreSpecies::AbsorbResultsFromWorkerScorer"<<G4endl;
193 {
194 SpeciesMap::iterator it_map1 = right->fSpeciesInfoPerTime.begin();
195 SpeciesMap::iterator end_map1 = right->fSpeciesInfoPerTime.end();
196
197 for (; it_map1 != end_map1; ++it_map1) {
198 InnerSpeciesMap& map2 = it_map1->second;
199 InnerSpeciesMap::iterator it_map2 = map2.begin();
200 InnerSpeciesMap::iterator end_map2 = map2.end();
201
202 for (; it_map2 != end_map2; ++it_map2) {
203 SpeciesInfo& molInfo = fSpeciesInfoPerTime[it_map1->first][it_map2->first];
204 molInfo.fNumber += it_map2->second.fNumber;
205 molInfo.fG += it_map2->second.fG;
206 molInfo.fG2 += it_map2->second.fG2;
207
208 // G4cout << "In AbsorbeResultsFromWorkerScorer: fNumber "
209 // << molInfo.fNumber
210 // << " fG "
211 // << molInfo.fG
212 // << G4endl;
213 }
214 }
215 }
216 //---------------------------------------------------------
217 #ifdef _ScoreSpecies_FOR_ALL_EVENTS
218 {
219 SpeciesMapPerEvent::iterator it_map1 = right->fSpeciesInfoPerEvent.begin();
220 SpeciesMapPerEvent::iterator end_map1 = right->fSpeciesInfoPerEvent.end();
221
222 for (; it_map1 != end_map1; ++it_map1) {
223 auto& map2 = it_map1->second;
224 InnerSpeciesMapPerEvent::iterator it_map2 = map2.begin();
225 InnerSpeciesMapPerEvent::iterator end_map2 = map2.end();
226
227 for (; it_map2 != end_map2; ++it_map2) {
228 SpeciesInfoSOA& molInfo = fSpeciesInfoPerEvent[it_map1->first][it_map2->first];
229 molInfo.fNumber.insert(molInfo.fNumber.end(), it_map2->second.fNumber.begin(),
230 it_map2->second.fNumber.end());
231 molInfo.fG.insert(molInfo.fG.end(), it_map2->second.fG.begin(), it_map2->second.fG.end());
232 molInfo.fG2.insert(molInfo.fG2.end(), it_map2->second.fG2.begin(),
233 it_map2->second.fG2.end());
234 molInfo.fEventID.insert(molInfo.fEventID.end(), it_map2->second.fEventID.begin(),
235 it_map2->second.fEventID.end());
236 // G4cout << "In AbsorbeResultsFromWorkerScorer: fNumber "
237 // << molInfo.fNumber
238 // << " fG "
239 // << molInfo.fG
240 // << G4endl;
241 }
242 }
243 right->fSpeciesInfoPerEvent.clear();
244 }
245 #endif
246
247 fNEvent += right->fNEvent;
248 right->fNEvent = 0;
249 right->fEdep = 0.;
250 }
251
252 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
253
254 void ScoreSpecies::DrawAll()
255 {
256 ;
257 }
258
259 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
260
261 void ScoreSpecies::PrintAll()
262 {
263 G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
264 G4cout << " PrimitiveScorer " << GetName() << G4endl;
265 G4cout << " Number of events " << fNEvent << G4endl;
266 G4cout << " Number of energy deposition recorded " << fEvtMap->entries() << G4endl;
267
268 for (auto itr : *fEvtMap->GetMap()) {
269 G4cout << " copy no.: " << itr.first << " energy deposit: " << *(itr.second) / GetUnitValue()
270 << " [" << GetUnit() << "]" << G4endl;
271 }
272 }
273
274 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
275
276 void ScoreSpecies::ASCII()
277 {
278 std::ofstream out("Species.Txt");
279 if (!out) return;
280
281 out << "Time is in ns" << G4endl;
282
283 for (auto it_map1 : fSpeciesInfoPerTime) {
284 InnerSpeciesMap& map2 = it_map1.second;
285
286 out << it_map1.first << G4endl;
287
288 for (auto it_map2 : map2) {
289 out << it_map2.first->GetName() << " " << it_map2.second.fNumber << G4endl;
290 }
291 }
292
293 out.close();
294 }
295
296 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
297
298 void ScoreSpecies::OutputAndClear()
299 {
300 if (G4Threading::IsWorkerThread()) return;
301
302 //----------------------------------------------------------------------------
303 // Save results
304
305 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
306 analysisManager->SetDefaultFileType(fOutputType);
307
308 if (analysisManager) {
309 this->WriteWithAnalysisManager(analysisManager);
310 }
311
312 fNEvent = 0;
313 fSpeciesInfoPerTime.clear();
314 }
315
316 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
317
318 void ScoreSpecies::WriteWithAnalysisManager(G4VAnalysisManager* analysisManager)
319 {
320 // G4cout << "ScoreSpecies::WriteWithAnalysisManager" << G4endl;
321 analysisManager->OpenFile("Species.root");
322 int fNtupleID = analysisManager->CreateNtuple("species", "species");
323 analysisManager->CreateNtupleIColumn(fNtupleID, "speciesID");
324 analysisManager->CreateNtupleIColumn(fNtupleID, "number");
325 analysisManager->CreateNtupleIColumn(fNtupleID, "nEvent");
326 analysisManager->CreateNtupleSColumn(fNtupleID, "speciesName");
327 analysisManager->CreateNtupleDColumn(fNtupleID, "time");
328 analysisManager->CreateNtupleDColumn(fNtupleID, "sumG");
329 analysisManager->CreateNtupleDColumn(fNtupleID, "sumG2");
330 analysisManager->FinishNtuple(fNtupleID);
331
332 for (auto it_map1 : fSpeciesInfoPerTime) {
333 InnerSpeciesMap& map2 = it_map1.second;
334
335 for (auto it_map2 : map2) {
336 double time = it_map1.first;
337 auto species = it_map2.first;
338 const G4String& name = species->GetName();
339 int molID = it_map2.first->GetMoleculeID();
340 int number = it_map2.second.fNumber;
341 double G = it_map2.second.fG;
342 double G2 = it_map2.second.fG2;
343
344 analysisManager->FillNtupleIColumn(fNtupleID, 0, molID); // MolID
345 analysisManager->FillNtupleIColumn(fNtupleID, 1, number); // Number
346 analysisManager->FillNtupleIColumn(fNtupleID, 2, fNEvent); // Total nb events
347 analysisManager->FillNtupleSColumn(fNtupleID, 3, name); // molName
348 analysisManager->FillNtupleDColumn(fNtupleID, 4, time); // time
349 analysisManager->FillNtupleDColumn(fNtupleID, 5, G); // G
350 analysisManager->FillNtupleDColumn(fNtupleID, 6, G2); // G2
351 analysisManager->AddNtupleRow(fNtupleID);
352 }
353 }
354
355 //----------------------------------------------------------------------------
356
357 #ifdef _ScoreSpecies_FOR_ALL_EVENTS
358 fNtupleID = analysisManager->CreateNtuple("species_all", "species_all");
359 analysisManager->CreateNtupleIColumn(fNtupleID, "speciesID");
360 analysisManager->CreateNtupleIColumn(fNtupleID, "number");
361 analysisManager->CreateNtupleIColumn(fNtupleID, "nEvent");
362 analysisManager->CreateNtupleSColumn(fNtupleID, "speciesName");
363 analysisManager->CreateNtupleDColumn(fNtupleID, "time");
364 analysisManager->CreateNtupleDColumn(fNtupleID, "G");
365 analysisManager->CreateNtupleDColumn(fNtupleID, "G2");
366 analysisManager->CreateNtupleIColumn(fNtupleID, "eventID");
367 analysisManager->FinishNtuple(fNtupleID);
368
369 for (auto it_map1 : fSpeciesInfoPerEvent) {
370 InnerSpeciesMapPerEvent& map2 = it_map1.second;
371
372 for (auto it_map2 : map2) {
373 double time = it_map1.first;
374 const Species& species = it_map2.first;
375 const G4String& name = species->GetName();
376 int molID = it_map2.first->GetMoleculeID();
377
378 size_t nG = it_map2.second.fG.size();
379
380 for (size_t i = 0; i < nG; ++i) {
381 int number = it_map2.second.fNumber[i];
382 double G = it_map2.second.fG[i];
383 double G2 = it_map2.second.fG2[i];
384 int eventID = it_map2.second.fEventID[i];
385
386 analysisManager->FillNtupleIColumn(fNtupleID, 0, molID); // MolID
387 analysisManager->FillNtupleIColumn(fNtupleID, 1, number); // Number
388 analysisManager->FillNtupleIColumn(fNtupleID, 2, fNEvent); // Total nb events
389 analysisManager->FillNtupleSColumn(fNtupleID, 3, name); // molName
390 analysisManager->FillNtupleDColumn(fNtupleID, 4, time); // time
391 analysisManager->FillNtupleDColumn(fNtupleID, 5, G); // G
392 analysisManager->FillNtupleDColumn(fNtupleID, 6, G2); // G2
393 analysisManager->FillNtupleIColumn(fNtupleID, 7, eventID); // EventID
394 analysisManager->AddNtupleRow(fNtupleID);
395 }
396 }
397 }
398 #endif
399
400 analysisManager->Write();
401 analysisManager->CloseFile();
402 }
403