BDSIM
BDSIM is a Geant4 extension toolkit for simulation of particle transport in accelerator beamlines.
BDSScorerMeshInfo.cc
1/*
2Beam Delivery Simulation (BDSIM) Copyright (C) Royal Holloway,
3University of London 2001 - 2022.
4
5This file is part of BDSIM.
6
7BDSIM is free software: you can redistribute it and/or modify
8it under the terms of the GNU General Public License as published
9by the Free Software Foundation version 3 of the License.
10
11BDSIM is distributed in the hope that it will be useful, but
12WITHOUT ANY WARRANTY; without even the implied warranty of
13MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14GNU General Public License for more details.
15
16You should have received a copy of the GNU General Public License
17along with BDSIM. If not, see <http://www.gnu.org/licenses/>.
18*/
19#include "BDSDebug.hh"
20#include "BDSException.hh"
21#include "BDSExtent.hh"
22#include "BDSScorerMeshInfo.hh"
23#include "BDSUtilities.hh"
24#ifdef USE_BOOST
25#include "BDSBH4DTypeDefs.hh"
26#endif
27
28#include "G4Types.hh"
29
30#include "parser/scorermesh.h"
31
32#include "CLHEP/Units/SystemOfUnits.h"
33
34#include <algorithm>
35#include <fstream>
36#include <iterator>
37#include <vector>
38
39BDSScorerMeshInfo::BDSScorerMeshInfo(const GMAD::ScorerMesh& mesh)
40{
41 name = G4String(mesh.name);
42 geometryType = BDS::LowerCase(G4String(mesh.geometryType));
43 nBinsX = mesh.nx;
44 nBinsY = mesh.ny;
45 nBinsZ = mesh.nz;
46 nBinsR = mesh.nr;
47 nBinsPhi = mesh.nphi;
48 nBinsE = mesh.ne;
49
50 if (geometryType == "box")
51 {
52 if (!BDS::IsFinite(mesh.xsize)) {
53 throw BDSException(__METHOD_NAME__, "xsize must be finite in mesh \"" + mesh.name + "\"");
54 }
55 if (!BDS::IsFinite(mesh.ysize)) {
56 throw BDSException(__METHOD_NAME__, "ysize must be finite in mesh \"" + mesh.name + "\"");
57 }
58 if (!BDS::IsFinite(mesh.zsize)) {
59 throw BDSException(__METHOD_NAME__, "zsize must be finite in mesh \"" + mesh.name + "\"");
60 }
61
62 if (!BDS::IsFinite(nBinsX)){
63 throw BDSException(__METHOD_NAME__, "nx must be finite in mesh \"" + mesh.name + "\"");
64 }
65 if (!BDS::IsFinite(nBinsY)){
66 throw BDSException(__METHOD_NAME__, "ny must be finite in mesh \"" + mesh.name + "\"");
67 }
68 if (!BDS::IsFinite(nBinsZ)){
69 throw BDSException(__METHOD_NAME__, "nz must be finite in mesh \"" + mesh.name + "\"");
70 }
71
72 }
73 else if (geometryType == "cylindrical")
74 {
75 if (!BDS::IsFinite(mesh.zsize)) {
76 throw BDSException(__METHOD_NAME__, "zsize must be finite in mesh \"" + mesh.name + "\"");
77 }
78 if (!BDS::IsFinite(mesh.rsize)) {
79 throw BDSException(__METHOD_NAME__, "rsize must be finite in mesh \"" + mesh.name + "\"");
80 }
81 if (!BDS::IsFinite(nBinsZ)){
82 throw BDSException(__METHOD_NAME__, "nz must be finite in mesh \"" + mesh.name + "\"");
83 }
84 if (!BDS::IsFinite(nBinsPhi)){
85 throw BDSException(__METHOD_NAME__, "nphi must be finite in mesh \"" + mesh.name + "\"");
86 }
87 if (!BDS::IsFinite(nBinsR)){
88 throw BDSException(__METHOD_NAME__, "nr must be finite in mesh \"" + mesh.name + "\"");
89 }
90
91 }
92
93 xLow = -0.5*mesh.xsize * CLHEP::m;
94 xHigh = 0.5*mesh.xsize * CLHEP::m;
95 yLow = -0.5*mesh.ysize * CLHEP::m;
96 yHigh = 0.5*mesh.ysize * CLHEP::m;
97 zLow = -0.5*mesh.zsize * CLHEP::m;
98 zHigh = 0.5*mesh.zsize * CLHEP::m;
99 rLow = 0 * CLHEP::m;
100 rHigh = mesh.rsize * CLHEP::m;
101 eLow = mesh.eLow* CLHEP::GeV;
102 eHigh = mesh.eHigh* CLHEP::GeV;
103 eScale = mesh.eScale;
104
105 extent = BDSExtent(xLow, xHigh,
106 yLow, yHigh,
107 zLow, zHigh);
108
109 if (eScale == "user")
110 {// In future we can move RBDS::BinLoader to a separate library and use that both here and in rebdsim
111 std::string const BinsEdgesFile(mesh.eBinsEdgesFilenamePath);
112 std::ifstream file(BinsEdgesFile.c_str());
113
114 if (file)
115 {
116 // Reading of the bins edges file.
117 std::istream_iterator<double> it(file);
118 std::istream_iterator<double> end;
119 std::back_insert_iterator<std::vector<double>> it2(eBinsEdges);
120
121 std::copy(it, end, it2);
122 }
123 else
124 {throw BDSException(__METHOD_NAME__, "eBinsEdgesFilenamePath must be the path to a .txt file");}
125
126 nBinsE = (G4int)eBinsEdges.size()-1;
127 eLow = eBinsEdges[0];
128 eHigh = eBinsEdges[nBinsE];
129 }
130
131 if (nBinsE > 1)
132 {
133#ifdef USE_BOOST
134 if (eScale == "linear")
135 {energyAxis = new boost_histogram_linear_axis(nBinsE, eLow, eHigh, "energy");}
136 else if (eScale == "log")
137 {energyAxis = new boost_histogram_log_axis(nBinsE, eLow, eHigh, "energy");}
138 else if (eScale == "user")
139 {
140 std::vector<double> eBinsEdgesEnergyAxis = eBinsEdges;
141 std::for_each(eBinsEdgesEnergyAxis.begin(), eBinsEdgesEnergyAxis.end(), [](double& el){el *= CLHEP::GeV;});
142 energyAxis = new boost_histogram_variable_axis(eBinsEdgesEnergyAxis, "energy");
143 }
144 else
145 {throw BDSException(__METHOD_NAME__, "eScale must be 'linear', 'log' or 'user' in mesh \"" + mesh.name + "\"");}
146#endif
147 }
148}
General exception with possible name of object and message.
Definition: BDSException.hh:35
Holder for +- extents in 3 dimensions.
Definition: BDSExtent.hh:39
ScorerMesh class for parser.
Definition: scorermesh.h:38
int ny
Number of bins in y.
Definition: scorermesh.h:45
int ne
Number of bins in E.
Definition: scorermesh.h:49
int nx
Number of bins in x.
Definition: scorermesh.h:44
int nphi
Number of bins in Phi.
Definition: scorermesh.h:48
int nz
Number of bins in z.
Definition: scorermesh.h:46
std::string eScale
E scaling type.
Definition: scorermesh.h:56
std::string name
Name of this placement.
Definition: scorermesh.h:40
double xsize
X total width.
Definition: scorermesh.h:50
double eHigh
E High limit.
Definition: scorermesh.h:55
int nr
Number of bins in R.
Definition: scorermesh.h:47
double eLow
E Low limit.
Definition: scorermesh.h:54
double ysize
Y total width.
Definition: scorermesh.h:51
std::string geometryType
Name of scorermesh geometry to use.
Definition: scorermesh.h:42
double rsize
R total length.
Definition: scorermesh.h:53
std::string eBinsEdgesFilenamePath
E bins edges filename path.
Definition: scorermesh.h:57
double zsize
Z total width.
Definition: scorermesh.h:52
G4String LowerCase(const G4String &str)
Utility function to simplify lots of syntax changes for pedantic g4 changes.
G4bool IsFinite(G4double value, G4double tolerance=std::numeric_limits< double >::epsilon())