| 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 /// \file scavenger/src/ScoreSpecies.cc
27 /// \brief Implementation of the scavenger::ScoreSpecies class
28
29 #include "ScoreSpecies.hh"
30
31 #include "G4AnalysisManager.hh"
32 #include "G4Event.hh"
33 #include "G4Scheduler.hh"
34 #include "G4UImessenger.hh"
35 #include "G4UnitsTable.hh"
36
37 #include <G4EventManager.hh>
38 #include <G4MolecularConfiguration.hh>
39 #include <G4MoleculeCounter.hh>
40 #include <G4SystemOfUnits.hh>
41 #include <globals.hh>
42
43 namespace scavenger
44 {
45
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
47
48 ScoreSpecies::ScoreSpecies(const G4String& name, const G4int& depth)
49 : G4VPrimitiveScorer(name, depth),
50 G4UImessenger(),
51 fpSpeciesdir(new G4UIdirectory("/scorer/species/")),
52 fpTimeBincmd(new G4UIcmdWithAnInteger("/scorer/species/nOfTimeBins", this)),
53 fpAddTimeToRecordcmd(new G4UIcmdWithADoubleAndUnit("/scorer/species/addTimeToRecord", this)),
54 fpSetResultsFileNameCmd(new G4UIcmdWithAString("/scorer/species/setRootFileName", this))
55 {
56 fpSpeciesdir->SetGuidance("ScoreSpecies commands");
57 fpSetResultsFileNameCmd->SetGuidance("Set the root file name");
58 G4MoleculeCounter::Instance()->ResetCounter();
59 }
60
61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
62
63 void ScoreSpecies::SetNewValue(G4UIcommand* command, G4String newValue)
64 {
65 if (command == fpAddTimeToRecordcmd.get()) {
66 G4double cmdTime = fpAddTimeToRecordcmd->GetNewDoubleValue(newValue);
67 AddTimeToRecord(cmdTime);
68 }
69 if (command == fpTimeBincmd.get()) {
70 ClearTimeToRecord();
71 G4int cmdBins = fpTimeBincmd->GetNewIntValue(newValue);
72 G4double timeMin = 1 * ps;
73 G4double timeMax = G4Scheduler::Instance()->GetEndTime() - timeMin;
74 G4double timeLogMin = std::log10(timeMin);
75 G4double timeLogMax = std::log10(timeMax);
76 for (G4int i = 0; i < cmdBins; i++) {
77 AddTimeToRecord(std::pow(10, timeLogMin + i * (timeLogMax - timeLogMin) / (cmdBins - 1)));
78 }
79 }
80 if (command == fpSetResultsFileNameCmd.get()) {
81 fRootFileName = newValue;
82 }
83 }
84
85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
86
87 G4bool ScoreSpecies::ProcessHits(G4Step* aStep, G4TouchableHistory*)
88 {
89 G4double edep = aStep->GetTotalEnergyDeposit();
90 if (edep == 0.) {
91 return FALSE;
92 }
93 edep *= aStep->GetPreStepPoint()->GetWeight(); // (Particle Weight)
94 auto index = GetIndex(aStep);
95 fEvtMap->add(index, edep);
96 fEdep += edep;
97 return TRUE;
98 }
99
100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
101
102 void ScoreSpecies::Initialize(G4HCofThisEvent* HCE)
103 {
104 fEvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName());
105 if (fHCID < 0) {
106 fHCID = GetCollectionID(0);
107 }
108 HCE->AddHitsCollection(fHCID, (G4VHitsCollection*)fEvtMap);
109 G4MoleculeCounter::Instance()->ResetCounter();
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 == nullptr || species->empty()) {
125 G4cout << "No molecule recorded, energy deposited= " << G4BestUnit(fEdep, "Energy") << G4endl;
126 ++fNEvent;
127 fEdep = 0.;
128 G4MoleculeCounter::Instance()->ResetCounter();
129 return;
130 }
131 for (auto molecule : *species) {
132 for (auto time_mol : fTimeToRecord) {
133 G4double n_mol = G4MoleculeCounter::Instance()->GetNMoleculesAtTime(molecule, time_mol);
134 if (n_mol < 0) {
135 G4cerr << "N molecules not valid < 0 " << G4endl;
136 G4Exception("", "N<0", FatalException, "");
137 }
138 SpeciesInfo& molInfo = fSpeciesInfoPerTime[time_mol][molecule];
139 molInfo.fNumber += n_mol;
140 G4double gValue = (n_mol / (fEdep / eV)) * 100.;
141 molInfo.fG += gValue;
142 molInfo.fG2 += gValue * gValue;
143 }
144 }
145 ++fNEvent;
146 fEdep = 0.;
147 G4MoleculeCounter::Instance()->ResetCounter();
148 }
149
150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
151
152 void ScoreSpecies::AbsorbResultsFromWorkerScorer(G4VPrimitiveScorer* workerScorer)
153 {
154 auto right = dynamic_cast<ScoreSpecies*>(dynamic_cast<G4VPrimitiveScorer*>(workerScorer));
155
156 if (right == nullptr) {
157 return;
158 }
159 if (right == this) {
160 return;
161 }
162
163 auto it_map1 = right->fSpeciesInfoPerTime.begin();
164 auto end_map1 = right->fSpeciesInfoPerTime.end();
165
166 for (; it_map1 != end_map1; ++it_map1) {
167 InnerSpeciesMap& map2 = it_map1->second;
168 auto it_map2 = map2.begin();
169 auto end_map2 = map2.end();
170
171 for (; it_map2 != end_map2; ++it_map2) {
172 SpeciesInfo& molInfo = fSpeciesInfoPerTime[it_map1->first][it_map2->first];
173 molInfo.fNumber += it_map2->second.fNumber;
174 molInfo.fG += it_map2->second.fG;
175 molInfo.fG2 += it_map2->second.fG2;
176 }
177 }
178 right->fSpeciesInfoPerTime.clear();
179 fNEvent += right->fNEvent;
180 right->fNEvent = 0;
181 right->fEdep = 0.;
182 }
183
184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
185
186 void ScoreSpecies::PrintAll()
187 {
188 G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
189 G4cout << " PrimitiveScorer " << GetName() << G4endl;
190 G4cout << " Number of events " << fNEvent << G4endl;
191 G4cout << " Number of energy deposition recorded " << fEvtMap->entries() << G4endl;
192
193 for (auto itr : *fEvtMap->GetMap()) {
194 G4cout << " copy no.: " << itr.first << " energy deposit: " << *(itr.second) / GetUnitValue()
195 << " [" << GetUnit() << "]" << G4endl;
196 }
197 }
198
199 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
200
201 void ScoreSpecies::OutputAndClear()
202 {
203 if (G4Threading::IsWorkerThread()) {
204 return;
205 }
206
207 auto analysisManager = G4AnalysisManager::Instance();
208 analysisManager->SetDefaultFileType(fOutputType);
209 this->WriteWithAnalysisManager(analysisManager);
210 fNEvent = 0;
211 fSpeciesInfoPerTime.clear();
212 }
213
214 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
215
216 void ScoreSpecies::WriteWithAnalysisManager(G4VAnalysisManager* analysisManager)
217 {
218 analysisManager->OpenFile(fRootFileName);
219 G4int fNtupleID = analysisManager->CreateNtuple("species", "species");
220 analysisManager->CreateNtupleIColumn(fNtupleID, "speciesID");
221 analysisManager->CreateNtupleIColumn(fNtupleID, "number");
222 analysisManager->CreateNtupleIColumn(fNtupleID, "nEvent");
223 analysisManager->CreateNtupleSColumn(fNtupleID, "speciesName");
224 analysisManager->CreateNtupleDColumn(fNtupleID, "time");
225 analysisManager->CreateNtupleDColumn(fNtupleID, "sumG");
226 analysisManager->CreateNtupleDColumn(fNtupleID, "sumG2");
227 analysisManager->FinishNtuple(fNtupleID);
228
229 for (const auto& it_map1 : fSpeciesInfoPerTime) {
230 const InnerSpeciesMap& map2 = it_map1.second;
231
232 for (const auto& it_map2 : map2) {
233 G4double time = it_map1.first;
234 auto species = it_map2.first;
235 const G4String& name = species->GetName();
236 auto molID = it_map2.first->GetMoleculeID();
237 auto number = it_map2.second.fNumber;
238 auto G = it_map2.second.fG;
239 auto G2 = it_map2.second.fG2;
240 analysisManager->FillNtupleIColumn(fNtupleID, 0, molID); // MolID
241 analysisManager->FillNtupleIColumn(fNtupleID, 1, number); // Number
242 analysisManager->FillNtupleIColumn(fNtupleID, 2, fNEvent); // Total nb events
243 analysisManager->FillNtupleSColumn(fNtupleID, 3, name); // molName
244 analysisManager->FillNtupleDColumn(fNtupleID, 4, time); // time
245 analysisManager->FillNtupleDColumn(fNtupleID, 5, G); // G
246 analysisManager->FillNtupleDColumn(fNtupleID, 6, G2); // G2
247 analysisManager->AddNtupleRow(fNtupleID);
248 }
249 }
250
251 analysisManager->Write();
252 analysisManager->CloseFile();
253 }
254
255 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
256
257 } // namespace scavenger