| 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 // Gorad (Geant4 Open-source Radiation Analysis and Design)
27 //
28 // Author : Makoto Asai (SLAC National Accelerator Laboratory)
29 //
30 // Development of Gorad is funded by NASA Johnson Space Center (JSC)
31 // under the contract NNJ15HK11B.
32 //
33 // ********************************************************************
34 //
35 // GRScoreWriter.hh
36 // Defines the printout format of primitive scorer
37 //
38 // History
39 // September 8th, 2020 : first implementation
40 //
41 // ********************************************************************
42
43 #include "GRScoreWriter.hh"
44
45 #include <map>
46 #include <fstream>
47
48 #include "G4MultiFunctionalDetector.hh"
49 #include "G4SDParticleFilter.hh"
50 #include "G4VPrimitiveScorer.hh"
51 #include "G4VScoringMesh.hh"
52
53 GRScoreWriter::GRScoreWriter()
54 {;}
55
56 GRScoreWriter::~GRScoreWriter()
57 {;}
58
59 void GRScoreWriter::DumpQuantityToFile(const G4String& psName,
60 const G4String& fileName,
61 const G4String& option) {
62
63 // change the option string into lowercase to the case-insensitive.
64 G4String opt = option;
65 std::transform(opt.begin(), opt.end(), opt.begin(), (int (*)(int))(tolower));
66
67 // confirm the option
68 if(opt.size() == 0) opt = "csv";
69 if(opt.find("csv") == std::string::npos &&
70 opt.find("sequence") == std::string::npos) {
71 G4cerr << "ERROR : DumpToFile : Unknown option -> "
72 << option << G4endl;
73 return;
74 }
75
76 // open the file
77 std::ofstream ofile(fileName);
78 if(!ofile) {
79 G4cerr << "ERROR : DumpToFile : File open error -> "
80 << fileName << G4endl;
81 return;
82 }
83 ofile << "# Mesh or volume name: " << fScoringMesh->GetWorldName() << G4endl;
84
85 using MeshScoreMap = G4VScoringMesh::MeshScoreMap;
86 // retrieve the map
87 MeshScoreMap fSMap = fScoringMesh->GetScoreMap();
88
89
90 MeshScoreMap::const_iterator msMapItr = fSMap.find(psName);
91 if(msMapItr == fSMap.end()) {
92 G4cerr << "ERROR : DumpToFile : Unknown quantity, \""
93 << psName << "\"." << G4endl;
94 return;
95 }
96
97 std::map<G4int, G4StatDouble*> * score = msMapItr->second->GetMap();
98 ofile << "# Primitive scorer name: " << msMapItr->first << G4endl;
99 if(fact!=1.0)
100 { ofile << "# Multiplication factor : " << fact << G4endl; }
101
102 G4double unitValue = fScoringMesh->GetPSUnitValue(psName);
103 G4String unit = fScoringMesh->GetPSUnit(psName);
104 G4String divisionAxisNames[3];
105 fScoringMesh->GetDivisionAxisNames(divisionAxisNames);
106 ofile << "# First three integer entries: index of a cell of a mesh, just one cell for a volume tally." << G4endl;
107 ofile << "# Forth entry: sum of scores." << G4endl;
108 ofile << "# Fifth entry: sum of squared scores." << G4endl;
109 ofile << "# Sixth entry: number of events with non-zero effect." << G4endl;
110 ofile << "# Seventh entry: relative statistical error in %." << G4endl << G4endl;
111 // index of the cell
112 ofile << "# i" << divisionAxisNames[0]
113 << ", i" << divisionAxisNames[1]
114 << ", i" << divisionAxisNames[2];
115 // unit of scored value
116 ofile << ", total(value) ";
117 if(unit.size() > 0) ofile << "[" << unit << "]";
118 ofile << ", total(value^2), number of entries, relative error (%)" << G4endl;
119
120 // "sequence" option: write header info
121 if(opt.find("sequence") != std::string::npos) {
122 ofile << fNMeshSegments[0] << " " << fNMeshSegments[1] << " " << fNMeshSegments[2]
123 << G4endl;
124 }
125
126 // write quantity
127 long count = 0;
128 ofile << std::setprecision(16); // for double value with 8 bytes
129 for(int x = 0; x < fNMeshSegments[0]; x++) {
130 for(int y = 0; y < fNMeshSegments[1]; y++) {
131 for(int z = 0; z < fNMeshSegments[2]; z++) {
132 G4int idx = GetIndex(x, y, z);
133
134 if(opt.find("csv") != std::string::npos)
135 ofile << x << "," << y << "," << z << ",";
136
137 std::map<G4int, G4StatDouble*>::iterator value = score->find(idx);
138 if(value == score->end()) {
139 ofile << 0. << "," << 0. << "," << 0;
140 } else {
141 G4double x1 = value->second->sum_wx()/unitValue*fact;
142 G4double x2 = value->second->sum_wx2()/unitValue/unitValue*fact*fact;
143 G4int n = value->second->n();
144 // rms is sigma = sqrt((x2-x1^2/n)/(n-1))
145 // var = mean +/- rms/sqrt(n): means that the relative error of the sum = rms*sqrt(n)/x1
146 G4double rms = value->second->rms()/unitValue*fact;
147 // Relative error in %
148 G4double relError = rms*std::sqrt(n)*100./x1;
149 ofile << x1 << ", " << x2 << ", " << n << ", " << relError;
150 ofile << G4endl;
151 G4double factor = relError*relError/100.;
152 if (factor > 1.)
153 {
154 G4cout << "# Mesh or volume name: " << fScoringMesh->GetWorldName()
155 << " -- # Primitive scorer name: " << msMapItr->first << G4endl
156 << " bin " << x << "," << y << "," << z << " : statistical error " << relError << "(%)" << G4endl
157 << " to reduce the statistical error below 10%, increase number of events approximately "
158 << factor << " times." << G4endl;
159 }
160 }
161
162 if(opt.find("csv") != std::string::npos) {
163 ofile << G4endl;
164 } else if(opt.find("sequence") != std::string::npos) {
165 ofile << " ";
166 if(count++%5 == 4) ofile << G4endl;
167 }
168
169 } // z
170 } // y
171 } // x
172 ofile << std::setprecision(6);
173
174 // close the file
175 ofile.close();
176 }
177
178 void GRScoreWriter::DumpAllQuantitiesToFile(const G4String& fileName,
179 const G4String& option) {
180
181 // change the option string into lowercase to the case-insensitive.
182 G4String opt = option;
183 std::transform(opt.begin(), opt.end(), opt.begin(), (int (*)(int))(tolower));
184
185 // confirm the option
186 if(opt.size() == 0) opt = "csv";
187 if(opt.find("csv") == std::string::npos &&
188 opt.find("sequence") == std::string::npos) {
189 G4cerr << "ERROR : DumpToFile : Unknown option -> "
190 << option << G4endl;
191 return;
192 }
193
194 // open the file
195 std::ofstream ofile(fileName);
196 if(!ofile) {
197 G4cerr << "ERROR : DumpToFile : File open error -> "
198 << fileName << G4endl;
199 return;
200 }
201 ofile << "# Mesh or volume name: " << fScoringMesh->GetWorldName() << G4endl;
202 if(fact!=1.0)
203 { ofile << "# Multiplication factor : " << fact << G4endl; }
204 ofile << "# First three integer entries: index of a cell of a mesh, just one cell for a volume tally." << G4endl;
205 ofile << "# Forth entry: sum of scores." << G4endl;
206 ofile << "# Fifth entry: sum of squared scores." << G4endl;
207 ofile << "# Sixth entry: number of events with non-zero effect." << G4endl;
208 ofile << "# Seventh entry: relative statistical error in %." << G4endl << G4endl;
209
210 // retrieve the map
211 using MeshScoreMap = G4VScoringMesh::MeshScoreMap;
212 MeshScoreMap fSMap = fScoringMesh->GetScoreMap();
213 MeshScoreMap::const_iterator msMapItr = fSMap.begin();
214 std::map<G4int, G4StatDouble*> * score;
215 for(; msMapItr != fSMap.end(); msMapItr++) {
216
217 G4String psname = msMapItr->first;
218
219 score = msMapItr->second->GetMap();
220 ofile << "# Primitive scorer name: " << msMapItr->first << G4endl;
221
222 G4double unitValue = fScoringMesh->GetPSUnitValue(psname);
223 G4String unit = fScoringMesh->GetPSUnit(psname);
224 G4String divisionAxisNames[3];
225 fScoringMesh->GetDivisionAxisNames(divisionAxisNames);
226 // index order
227 ofile << "# i" << divisionAxisNames[0]
228 << ", i" << divisionAxisNames[1]
229 << ", i" << divisionAxisNames[2];
230 // unit of scored value
231 ofile << ", total(value) ";
232 if(unit.size() > 0) ofile << "[" << unit << "]";
233 ofile << ", total(value^2), number of entries, relative error (%)" << G4endl;
234
235
236 // "sequence" option: write header info
237 if(opt.find("sequence") != std::string::npos) {
238 ofile << fNMeshSegments[0] << " " << fNMeshSegments[1] << " " << fNMeshSegments[2]
239 << G4endl;
240 }
241
242 // write quantity
243 long count = 0;
244 ofile << std::setprecision(16); // for double value with 8 bytes
245 for(int x = 0; x < fNMeshSegments[0]; x++) {
246 for(int y = 0; y < fNMeshSegments[1]; y++) {
247 for(int z = 0; z < fNMeshSegments[2]; z++) {
248 G4int idx = GetIndex(x, y, z);
249
250 if(opt.find("csv") != std::string::npos)
251 ofile << x << "," << y << "," << z << ",";
252
253 std::map<G4int, G4StatDouble*>::iterator value = score->find(idx);
254 if(value == score->end()) {
255 ofile << 0. << "," << 0. << "," << 0;
256 } else {
257 G4double x1 = value->second->sum_wx()/unitValue*fact;
258 G4double x2 = value->second->sum_wx2()/unitValue/unitValue*fact*fact;
259 G4int n = value->second->n();
260 // rms is sigma = sqrt((x2-x1^2/n)/(n-1))
261 // var = mean +/- rms/sqrt(n): means that the relative error of the sum = rms*sqrt(n)/x1
262 G4double rms = value->second->rms()/unitValue*fact;
263 // Relative error in %
264 G4double relError = rms*std::sqrt(n)*100./x1;
265 ofile << x1 << ", " << x2 << ", " << n << ", " << relError;
266 ofile << G4endl;
267 G4double factor = relError*relError/100.;
268 if (factor > 1.)
269 {
270 G4cout << "# Mesh or volume name: " << fScoringMesh->GetWorldName()
271 << " -- # Primitive scorer name: " << msMapItr->first << G4endl
272 << " bin " << x << "," << y << "," << z << " : statistical error " << relError << "(%)" << G4endl
273 << " to reduce the statistical error below 10%, increase number of events approximately "
274 << factor << " times." << G4endl;
275 }
276 }
277
278 if(opt.find("csv") != std::string::npos) {
279 ofile << G4endl;
280 } else if(opt.find("sequence") != std::string::npos) {
281 ofile << " ";
282 if(count++%5 == 4) ofile << G4endl;
283 }
284
285 } // z
286 } // y
287 } // x
288 ofile << std::setprecision(6);
289
290 } // for(; msMapItr ....)
291
292 // close the file
293 ofile.close();
294
295 }
296
297
298